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On the Numerical Solution of Sturm-Liouville 
Differential Equations 


1, Introduction. It is the purpose of this paper to use the well-known relation 
that exists between a Sturm-Liouville differential equation together with its 
boundary conditions and normalization condition and a problem in the calculus 
of variations, to obtain an algebraic characteristic value problem from whose 
solution we may obtain approximations to the characteristic values and character- 
istic functions of the Sturm-Liouville equation. That is, we shall obtain approxi- 
mations to the values of \ and the functions y(x) satisfying the equation 


(1.1) y” + (q(x) + Ar(x))y = 0, 
such that the boundary conditions 


ay (a) + asy’(a) = 0 


(1.2) 
Bry(b) + Bxy’(b) = 0 


and the normalization condition 


(1.3) f ‘ rydx =1 


are satisfied, where the prime denotes the derivative of y(x) with respect to x, 
q and r being functions of x, and the latter function is required to be positive 
over the interval 

a<x<b. 


Equations of the form 


(py’)’ + gy + Ary = 0 


may be reduced to that of the form of equation (1.1) by a transformation of the 
independent variable. 

The equations (1.1) and (1.2) themselves may be approximated by a set of 
linear algebraic equations as follows (cf., Hildebrand [1]). Let us use a set of 
discrete values of x 


x =atth (¢ =0,1,---,N +1) 


b-a 


o "Fas 


and define 
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Then we may approximate y’(x;) = y,’ and y’’(x;) = y,/’ by means of finite differ- 
ences of the y;. We would then obtain a set of equations of the form 


(1.4) Lis = 19: 


for the unknowns J and y;. That is, we would thus reduce the system of equations 
(1.1) and (1.2) to an algebraic characteristic value problem. However the matrix 
L,; would not in general be symmetric. Indeed as the order of the differences used 
to replace the derivatives increases, the asymmetry of the matrix L;; increases. 
This asymmetry occurs because the same type of difference formula cannot be 
used at the ends of the interval (a,b) as in the interior of the interval. 

If the matrix L;; is not symmetric, the reality of the values of \ cannot be 
guaranteed. Moreover the determination of the values of \ such that equation 
(1.4) has a solution is a much greater computational problem when L,; is an 
asymmetric matrix than when it is a symmetric one. The methods that will be 
used in the sequel to replace the system of equations (1.1) to (1.3) by an algebraic 
characteristic value problem will always lead to symmetric matrices. 

2. The Variational Problem. It is well known that equation (1.1) is the Euler 
condition that the integral 


(2.1) f ” (y)? — aye 


be an extremal where the functions y(x) satisfy equation (1.3). 
Hence the system of equations (1.1) to (1.3) is equivalent to finding functions 
y(x) subject to equations (1.2) such that the integral 


(2.2) 10) = [Co — aide - a [ras 


is an extremum where J is a Lagrange multiplier and the functions y(x) are re- 
stricted to satisfy the conditions (1.2). 

Instead of dealing with the differential equation (1.1) we deal with the integral 
(2.2). There are two methods that we shall consider: (1) approximating this 
integral by using quadrature and finite difference formulas (Method I), and (2) 
restricting the class of admissible functions to be made up of polynomial arcs 
(Method II). The latter method is of course a variant of the Rayleigh-Ritz 
method. However it differs from that method in that it uses a set of polynomials, 
each of degree m, say, and each one defined over an interval x; < x < Xi41, to 
represent an admissible function y(x) instead of an expansion of y(x) in terms of 
a convenient orthonormal set of functions. The latter representation is required 
to hold over the entire interval (a,b) and many terms may be needed because of 
the behavior of y(x) on a few of the intervals x; < x < xj41; whereas in the former 
representation we may use intervals of varying sizes or allow the admissible 
function to be made up of polynomial arcs, where the degree of each arc need 
not be the same. 

3. Method I. In the discussion of this method we shall restrict ourselves to 
the case where a2 = 82 = 0 in equations (1.2). The discussion may be readily 
modified to take care of more general boundary conditions. 
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lifter. The first step is to obtain a finite difference approximation to (2.2). To do 
this we introduce the quadrature formula 
b N+1 
tions (3.1) f o(x)dx = » Wid; + €1, 
atrix ty 
used where ¢; is the error in this formula. The weighting factors w; for the functional 
ases. values ¢; are to be determined later. We also introduce the formula 
ot be 
N41 
(3.2) y! = Ir ya Oifys + 15 ( = 0,1, --+,-N +1) 
ot be _— 
ation 
is an which is a finite difference approximation to the first derivative at x; in terms of 
ill be the functional values y;. The values of a;; are to be determined later also. 
braic Using (3.1) and (3.2) we approximate the first integral in (2.2) obtaining 
euler ree a“ re . 
L(y’)? — gy Idx = DY willy’)? — qa?#] +a 
a i=0 
N+1 N+1 N+1 
= 2 wih" LD agjt+nP -— LD way? +a 
i=0 7=0 i=0 
° N+1 N+1 
Gons = 2 [*) waits — waa lim + (4 + @) 
i,k=O i=0 
where 
N+1 N+1 
ee = DL wl2nh" DL auy; + 27] 
e re- i=0 7=0 
' and 6, is the Kronecker delta. Our equation becomes 
oral 
this » N+1 
| (2) 3) f L(y’)? — gy idx = SY Livin + (2 + &) 
arcs . — é; 
Ritz if we set 
ials, 
N+1 
1, to (3.4) Liz = > Wii — W3Q) je- 
is of i<0 
\ired é 
e of Likewise the second integral in (2.2) becomes 
mer ; pi? 
sible (3.5) f rydx = ys W5? Dj iVe + €3. 
need a i,k=0 
_ Thus I(y) is approximated by the quadratic form 
dily N41 N+1 


(3.6) F= > Liye —d* XL wrdayin, 


i,k=0 7,k=0 
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and if we insert the boundary conditions yo = yw+1 = 0 we obtain simply 


N N 
F= > Lpyim —d* LD wyrdpayinr. 
j,k=1 j,k=1 


The second step in our numerical method is to minimize the quadratic form. 
This is a simple problem in the calculus and the necessary condition for a minimum 
is that 


. ee = 

—— = L Lay — * LD wrdam = 9, 

OY; k=1 k=l 
for 7 = 1,---,N. This system of linear equations can be written in matrix 
notation 
(3.7) Ly — »\*WRy = 0, 


where L is a symmetric N X N matrix and WR is a diagonal one. 
In order to put (3.7) in a form analogous to (1.1) we multiply through by 
—1/w;, and if we define tj, to be 


tu = — (1/w,)Li, 
we can write 


N 
LX tive + A*ry;5 = 0, 


k=1 
which becomes 
Ty + A*Ry = 0, 


in matrix form. Using (3.4) and setting 


G; = 1/w;h?, 
we can write i, as 
N+1 


(3.8) te = — Gi D waits + gdp. 
i—0 


We note that Z and 7 are functions of w; and a,; from the approximation 
formulas (3.1) and (3.2) respectively. In the next section a method will be 
discussed for selecting w; and a;; so as to have T approximate the differential 
operator . 


a? 
D =—-+ ’ 
aa * 


to a specified order of approximation. 

4. The Order of Approximation of the Method. It was stated in the last 
section that the matrix operator T is to be an approximation to the differential 
operator D,. We would like to specify the order of this approximation. However, 
before doing so let us define what we mean by a certain “order of approximation” 
in this case. 
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In order to compare a differential operator with a matrix operator we need 
to arrive at a common basis for comparison. To do this consider the vector 


| qyityn” , 


m1 =] gostyi” 








| gnyntyn” 


whose components are obtained by evaluating Dzy at the discrete set of mesh 
points x;, 7 = 1, ---, N. Also consider the vector 
| N 
DX ture 


k=1 


N 
ve =| > tare 
k=1 








N 
X tweye 
ro 
Definition. If 
N 

~y X tan = ays + 94" + 00H") 
for 7 = 1, ---, N, then the matrix T approximates the differential operator D: 
to order o. 


This definition can be used as a basis for the selection of w; and a,; for formulas 
(3.1) and (3.2) respectively. To see this let us require that w; and a,; be chosen in 
such a way that the numerical method is of order c. Let us further require that 


N+1 


(4.2) yi! =h" YS aan + 0(h"), 
k=0 


ie., that (3.2) be an approximation formula also of order ¢. This means that 
the a, must satisfy the equations 


N+1 


(4.3) > @a(k —1)*=6, (wu = 0,1, ---,¢). 
k=0 
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For, if we expand y; about y;, we get 
@ ; he 
m= DL (k-ty™— 
u=0 B- 


and thus, in (4.2), 





N+1 © f N+ pr 
kD aam = LD | DL au(k — iy] — 
k=0 u=0 k=0 M- 
from which (4.3) follows. 
Using the same procedure, and recalling that yo = yyw41 = 0, we can write 
the left side of equation (4.1) as 


N N+1 N+1 N+1 
LD tame = DL tere = {y,[ ta| + ma | Dd (k - ta 
k=1 0 k=0 k=0 


h? N+1 
toi" | x (k ~ tn | + vf, 


k=0 


Thus we have from equation (3.8) 


N N+1 N+1 
YL tan = ys | LD (-G; LD waigu + asa)| 
k=0 i=0 


k=1 


N41 N41 
+ yi] LD (kk -j)(-G; XL waia + aa) | 
k=0 - i=0 


yep N+ w+ 
+o | xo kk -j(-G XL wan + ass)| 
i=0 


2!L imo 


Now we are in a position to treat (4.1) as an identity~in y; and its derivatives 
by using the last equation. By equating coefficients of y; and its derivatives up 
to y°t), we obtain 


N+1 N+1 


XL (-— GL wags + qin) = 45, 
k=0 =0 
N+1 N+1 if 


LX kk -7A(-—G; LX wagon + gon) = 0, 
k=0 i=) 


2 N+1 N+1 


~>ZXL (k-jP(-—G; XL wan + gon) = 1, 
2 k=0 i=0 


and 


N+1 N+1 


LX (k- A (-G X waa + gn) = 0, 
k=0 i=0 


; 
e 


wig 


} 
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: where tr = 3, ---,o + 1. Since 
5 N+1 
E LX qin = Gi, 
ie k=0 
r 
and 
i (k ‘i Doin = 0, 
and using the definition of G;, these become 
| N+1 N+1 
: p W Ai; > axn=0, 
ite i=0 k=O 
) N+1 N+1 
ya Wi; Zz, an(k — j) = 0, 
i=0 k=0 
j N+1 N+1 
> wai; + aun(k — 7)? = —2w,;, 
, P 1=0 k=0 
and 
N+1 N+1 
LX way L au(k — 7)" = 9, 
7=0 k=0 


where tr = 3, ---,o +1. 

A further simplification results if we use equation (4.3), i.e., if we make use 
of the fact that the a;; give us an approximation to the derivative of order o. 
For we can write 


»| ¥ aalk —j" = ¥ aal(k-J + GAP 





N+1 


D aal(k — i + u(k — GE —§) +--+ alk - DG - jet G- HI 


ives 

up 
= p(t — jr. 

With this result our equations become 


N+1 
) wai = 0, 
i=0 


N+1 


LY waiys(j — 1) 
i=0 


W;, 


and 


N+1 
¥ wa;(j —i)? =0 (7 =2,---,0), 
=0 
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or simply 
N+1 

(4.4) > w0i;(j —1)?’ =wrr (7r = 0,---, 0). 
t=0 

In these equations j = 1, ---, NV. 


Thus if we select a,;; and w; satisfying equations (4.3) and (4.4) our method 
will give an approximation of order o. 

5. The First Order Approximation. In this case ¢ = 1 and equations (4.3) 
and (4.4) become 


N+1 


LX ae = 0, 
(5.1) k=0 
N+1 
> an(k —7) =1 (@ =0---N +1), 
k=0 
and 
N+1 
bm Wii = 0, 
(5.2) imo 
N+1 
=w (j=1,---,N), 


LY wai — 1) 


respectively. Note that equations (5.1) give conditions on the rows of the matrix 
(ax) whereas equations (5.2) give conditions on the columns (except for the first 
and last) of the same matrix. 

It is easily verified that the matrix 








—1 1 0 0 0 
0 -1 1 0 0 
0 0 -1 1 0 
0 0 0 -1 1 . 
(5.3) (aa) = " : ; ; ; ; ‘ : : hy 
-_ —1 1 0 0 
0 -1 1 0 
, 0 0 -1 1 
I pt mgpinigg pg uesyg 
and w; satisfying ' 
(5.4) ~=s 


Wi = Wi41 (¢ =0,---,N —1) 


constitute a solution of (5.1) and (5.2), ie., equations (5.3) and (5.4) give us 
a and w; for the first order case. 

Equation (5.3) states that we use the simple first order forward difference 
formula 


yl = (viens — ¥s)/he + Lyi (h/2!) + 95 (2/3!) + +--+], 
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for i = 0, ---, N, and the simple first order backward difference formula 
ye = (yi = yi-1)/h = Lyi" (h/2!) _ yi (h?/3!) + ree], 


fori = N + 1. Equations (5.4) indicate that the quadrature formula to use is the 
simple first order ‘‘Rectangle Rule” 


. N 
(5.5) f f(x)dx = wi & fi — [fr'(h/2) + ---). 
a i=0 


Obviously 
w, = (6 — a)/(N + 1) = 4, 


since (5.5) must hold for f(x) = 1. 
6. The Second Order Approximation. In this case ¢ = 2 and equations (4.3) 
and (4.4) become 


N+1 


7. tx = 0, 
k=0 


N+1 


(6.1) LD aa(k — i) = 1, 
k=0 


N+1 
E aalk—1)? =0 ( =0,---,N +1), 
k=0 
and 
N+1 
LX wWaij = 0, 
1=0 


N+1 


(6.2) LD wais(j —1 


=0 


~~ 
— 
II 


Wj, 


N+1 


L way(G-1?=0 (G=1,---,), 
i=0 


respectively. 
It has been shown by Gregory [2] that when N is odd a solution to equations 
(6.1) and (6.2) is given by the matrix (aa) of the form 


[_—3/2 2-1/2 0 0 0 0 1 
-1/2 0 1/2 0 0O 0 0 
1/4 -1 0 1-1/4 0 0 
0 0-1/2 0 1/2 O 0 
scEncrcter Me aetatie-serk abr oid ieee aa 
0 -1/2 0 1/2 0 0 
0 1/4 -1 0 1-1/4 
0 0 0 --1/2 0 1/2 
; 0 0 0 1/2 -2 3/2] 
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and the w; are given by 


(6.4) 


WwW = 4wp, 
We. => 2wo, 
Ws = 4wo, 


wri = 2wWo, 
Wy = 4wWo, 
Wn+1 = Wo. 





If N is even the solutions differ only in the last four rows of the matrix (a) 
and the last four values of w,, i.e., 


(6.5) @)= 


and 


(6.6) 





-$/2 2-82 

~1/2 0 

1/4 -1 0 
0 
0 
0 
0 
0 


1/2 


0 0 0 
0 0 0 
1 —1/4 0 
—1/2 0 1/2 0 0 
4/17 —16/17 —3/34 18/17 —9/34 
0 0 —1/3 -—1/2 1 
0 ® -14 oO =1/2 
0 0 1/3 -1/2 -1 
W, = 4w», 
We = 2wWo, ° 


Wn-4 = 2wWo, 
Wn-3 = 4Wo, 


17 
Wn-2 = 7 Wo, 
27 
wWy-1 = rs Wo, 
27 
wn = ry Wo, 
9 
Wn+i1 = 3 


Oo |, 
0 
—1/6 
2/3 

7/6 | 





Equations (6.6) can be interpreted as a quadrature formula consisting of Simp- 
son’s Rule for all but the last three intervals and the ‘3/8 Rule’”’ for the last three 
intervals. This special treatment at the end is due to the fact that Simpson’s 
Rule can only be used with an odd number of mesh points and here JN is even. 

In a later section the solutions of equations (3.7), where L;; is defined by 
equation (3.4) for the first and second order approximations, will be compared 
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to the characteristic values and values of the characteristic functions at the mesh 
points of the Sturm-Liouville system described by the equation 


dy 
(6.7) a8 +rAy = 0 
with the boundary conditions 


(6.8) y(0) = y(1) = 0 


and the normalization condition 


(6.9) f yt dx = 1, 


This simple example for which the solutions are known to be 
(6.10) An = 2°n? 
(6.11) yn = ¥2 sin xnx 


will be used in the discussion of the effectiveness of Methods I and II. 

7. Method II. We now turn our attention to the derivation of the algebraic 
characteristic value problem associated with Method II. In this method, instead 
of approximating the quadratures occurring in (2.2) by mechanical formulas, the 
unknown solution (x) itself is replaced by an approximating function Y(x). The 
approximating function Y(x) is defined in terms of N + 2 parameters y; in a 
manner such that Y(x) assumes the value y; at the mesh point x;, i.e., Y(x;) = y; 
(j = 0, ---, N+ 1). More specifically, in each interval x; < x < x;4; the ap- 
proximating function Y (x) is defined to be a polynomial Y;(x) of degree m passing 
through y; and m other prescribed nearby ordinates. The boundary conditions 
(1.2) determine two of the ordinates in terms of the others. Then the integral 
(cf., equation (2.2)) 


(7.1) 1(Y) = 4: (¥"? — g¥? — rr ¥*)dx 


is a function of the values y; assumed by Y(x) at N distinct selected points. Re- 
quiring that the approximating function Y(x) be an extremal of (7.1) leads, 
without further approximations, to an algebraic characteristic value problem de- 
termining N sets of values for \ and the ordinates y;. 

The polynomials Y;(x) defining the function Y(x) are given by the Lagrange 
interpolation formula 


a2) Ye) 2 aso @— 80) 


(xi — %4) -** (eq — in) 





Vig + coe 


(w — x4) --- (@—%,) 
(xi, — %%) °°* (Xin — Zina)” 





+ 
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where % = 7+ k — 1; (k = 0,1, ---, m) and the /; are chosen such that Y(x) is 
given by “forward” or ‘“‘backward”’ interpolation formulas near the end points 
and by nearly ‘“‘central’’ formulas at interior points and such that the interpolation 
polynomials are chosen symmetrically with respect to the midpoint (a + 8)/2. 
More specifically, in obtaining the numerical results reported in this paper NV 
is odd and the choice of the ordinates through which each Y;(x) is required to 
pass is equivalent to the following rules determining the single integer /; 


l= fori < 34m -—1 
(7.3) 3m —2<1;< 4m—1 for }m—1<i1<4N 
2 am <1; <4m+1 for 3N <i<N+1-— im 
l=it+tm-N-1 fori > N +1 — 3m. 


Henceforth we shall assume that the mesh points x; are equally spaced with 
Xo = @, Xw41 = band x41 — x; = h = (6 — a)/(N +1). 

Since each Y;(x) depends linearly upon the ordinates y,, its derivatives do 
also. Hence we may write for the rth derivative of Y;(x) at x; 


hr ld 
(7.4) r! Yi (x) = L ati. 


j=l 


For fixed values of 2 and r, at most m + 1 of the coefficients a',; are different 
from zero. The summation in (7.4) extends from 7 = 1 to 7 = N since yo and 
yn+1 may be determined in terms of the others by the boundary conditions (1.2). 
The a‘,; thus will in general depend upon the boundary conditions. In the special 
case where az = 82 = 0, however, we have yo = yi: = 0, and the a‘,; are rational 
numbers which can be expressed as integers divided by m! From (7.4) it follows 
by Taylor’s theorem that 


(7.5) VG) EP ate - 2y. 


r=0 j=1 


We shall discuss the determination of the quantities a‘,; in the next section. 
In the remainder of this section we assume that they are known and evaluate 
equation (7.1) for 


(7.6) Y(x) = Vilx) x Sx < x41 m 


where Y;(x) is given by equation (7.5). 
Substituting equation (7.6) into (7.1) gives 


1(¥) = Ef" (LEE rhratyle — 2} 
— [La(e) + r(#) ILE EX bratsgys x — x5)" Pde. 


(7.7) 





a’; 


con 
titi 


Sin 
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is If we define 
its a ws 
on (7.8) Pi, = mah (x — x,)"+** dx, 
y 3 
N 1 Zzi+1 
to (7.9) Oi, = ova J. q(x) (x — x,)"+* dx, 
; 1 zi+1 
(7.10) Ri, = an f r(x) (x — x,)'+ dx, 
Nom 
(7.11) Li = 7 > ai, (P+, _ Ot )a's;, 
i=0 r,t=0 
th N m 
(7.12) Da = > aR in, 
do i—O 1, tm 


equation (7.7) becomes 


N 
I(Y) = DL Liz — ADix) yi. 
j,k=0 


The extrema of J(Y) are determined by the numbers y; satisfying 


nt 
nd N N 
2). (7.13) LD Lage =D Dye. 
ial k=l k=l 
l , . , . : ; 
Pe for 7 = 1, ---, N. This equation may be written in matrix notation as 
(7.14) LY = »DY 
when Y is the matrix of a single column: 
¥1 
mn. Y = 
ite 
YN 
8. Determination of a‘,;, Formulas for the determination of the numbers 
a',; may be obtained by differentiation of equation (7.2) and using the boundary 
conditions. However, as will be shown below we may also determine these quan- 
tities by a double induction, one on the index i and one on the degree of the 
polynomials involved. For simplicity in the ensuing discussion we assume that 
we are dealing with the special boundary conditions 
(8.1) Fo *@ Feu = 0. 
1x. 


Similar methods apply in the case of general boundary conditions. 
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We first discuss the inductive procedure we use to determine the a’; which 
are defined by the equation 


hrVo (xo) Xt} 

(8.2) a =) oi 

r! pars 
where Yo" (xo) is the rth derivative of the polynomial Yo(x) evaluated at x = xp. 
By definition Yo(x) is the mth degree Lagrange polynomial passing through 
Yo, V1, bec” Vm- 

Let P,(x) be the sth degree Lagrange polynomial passing through yo, yi, -- -, ¥.. 
Then 

Yo(x) = Pa(x). 

We note that 


(8.3) P(x) = Vo. 


In general we have 





(8.4) P,(x) = LD Par (x)% 
k=0 
where 
_  (*& — Xo) --- & — Xn-1)(% — Xeq1) °° (K — Xe) ; 
(8.5) Pale) = (x~ — X0) +++ (Xe — Xe—1) (Xe — Hep) -* > (Xe — Xs) 


We may also write 


“ PP, 0 
P(x) = DY 


r=0 


(x — Xo)" 


where P, ” (xo) is the rth derivative of P,(x) evaluated at x = xo. We may further 
write 


hr N+1 
(8.6) — P. (x0) = DL Beye 
r. k=0 
Then 
(8.7) a; _ B™,; 
and 
By = = 09 (a 
(8.8) aa Psy (Xo). 


It follows from equations (8.4) and (8.5) that 


(8.9) Pers = ——"" Pa b<stl 


Xk —~ Xe+1 
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fe (x — Xs) (x. — x0) (%_ — %1) --- (X_ — Xe-1) 
Pee, Pen et GO al Gn = eee jena 
_ & — *) 
he +1) a7 4) Pool), 


the last equation holding when the intervals are equal in length. 
Differentiating equations (8.9) and (8.10) r times we have 











Le ’ = 
PS 414 (x) = peeps b (x) w= = Px” (x) 
Xs41 
(r) ih (7) _ prt) 
Ps41s4i (x) = roe rs c ss (X) + hs . 1) Pss_ (x). 
In view of equation (8.8) these equations may be written as 
+1 1 
At ‘tf = — k= = 
(8.11) hit = he - 
s 1 
oie = — ——- fh + ——- Ba > 


s+1°” s+i1 
It follows from equation (8.3) that 


Boo = 1 
a =0 rX#0,k ¥0. 


These equations together with equations (8.11) define 6; and hence ay in view 
of equation (8.7). 

To obtain formulas by which the a}; may be computed for i > 0 we consider 
first the case in which Y;,(x) and Yi4:(x) are the same polynomial. The aj; and 
the aj. i. then determine the derivatives of this polynomial at x; and xj: respec- 
tively. Since Y;(x) = Yiy1(x) we have in particular 


r 


Fu. h’ 4 
a VY!” (x41) = of VY (xi41). 





Now 
ht 
1") (4. - (t) 
Y; (X41) r Y; ) GE —r)! 
hr’ » t! ht 
aie a pears... same... () 
3 VY}? (xi41) = 2 HG it Y;" (x;) 
m t! N 
- ~ Leni L ais. 


h’ —~ar 
Since = _ (Wh = x aj;*'y; it follows that 
7=1 
ri(t — rv)! 


(8.12) ah 


t=r 
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If N is chosen to be an odd integer and if the ordinates through which each Y;(x) 
is required to pass are chosen in accordance with (7.3) then (8.12) is applicable 
for 1 < $m — 2 and fori > N +1 — 3m. 

When the ordinates through which Y;(x) and Y;,:(x) are chosen in the same 
way, as for example when Y;,:(x) passes through y;:-1, yi, and yiz1 and Y;,;(x) 
passes through 4;, yi+1, and yi+2, then clearly 

anf = af 5-4. 
This formula is applicable, in accordance with (7.3), except near the end points 
and at the midpoint xy, where M = 3(N + 1). If m is even we find that Yy(x) 
= Yy-:(x) so that (8.12) applies. If m is odd then it follows from (7.3) that 
Yu(x) = Yu_—2(x) so that the induction formula (8.12) must be used twice. Noting 
that aMZi = a/f~? we have 
m s! m t! 


M n ° ~ 


ri(is—r)lislt — 





s=r t=s 


when m is odd. Summarizing these results 





m 
ta 2 
a? 
a for i>N+1—-— 
° = eee oe —_—— 
~~ = ri(t — vr)! = 2 
1 = M —1, if m is even 
“= eee Tee 
_ = : 1 1S 
 ” Sis — = 1¢ -— . 
m™ 
37 2<t<M-1 
at) = ab 5-1 for 


M-2<i<N+1-2 


where, again, N is odd and M = }(N + 1). These equations together with our 
previous results for determining the a°,; constitute a means for computing the 
coefficients a‘,; forz = 0, ---, N;r = 0, ---,m; andj = 1, ---, N for the bound- 
ary conditions y(a) = y(b) = 0 when N is chosen to be odd. 

9. The Computation Procedure. Both Method I and Method II lead to alge- 
braic problems of the type given by equation (7.13). In the former case the 
matrix D is diagonal and various simplifications occur. For both methods the 
computation proceeds in two similar steps: (1) computing the elements of the 
matrices L and D, and (2) solving equations (7.13) for \ and the y;. 

In Method I the matrix L is determined from the matrices a;; given by equa- 
tions (5.3) or (6.3) in accordance with equation (3.4). 

In Method II step (1) begins with the computation of a®,;. We have initially 
a matrix with m + 1 columns, r = 0, ---,m and N +1 rows, j = 0,---,N, 
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only the element (0,0) being different from 0. Then the transpose of the matrix 
\|a°,;|| is computed using the formulas derived in section 8. Thereafter the first 
row is no longer needed (assuming boundary conditions y(a) = y(b) = 0). Then 
we have routines which perform the following functions: (1) compute |{a‘,,;|| from 
\lai;"*|| using equations (8.13) ; (2) compute the matrix ||P‘ — Q*|| using equations 
(7.8) and (7.9); (3) matrix multiplication So a‘:(P‘.— Q+.); (4) matrix 
multiplication }70 a@‘x[_>,70 arj(P ‘+: — Q*++)] and addition of the result to the 
previous partial sum; (5) compute the matrix ||R‘|| using equations (7.10); 
(6) (master routine) step 7 and repeat until = N + 1. Routines (3) and (4) are 
used to form the partial sums for both L and D. Likewise (3) and (5) may con- 
veniently be combined into a single routine. Since || P* — Ql], ||R*|| are both sym- 
metric and can use the same storage locations, since L and D are symmetric, and 
since |\a‘||7-||P* — Q‘|| and |la*‘||7-||R*|| can use the same storage locations,-a 
total of 3(m + 1)(m + 2) + N(N + 1) + 2N(m + 1) matrix storage locations 
are required. It is also time-saving to store the integers 0! through m! for use in 
computing the ||a‘l|. 

If the functions g(x) and r(x) are symmetric with respect to x = 4(a + 6) 
it is time-saving to choose N to be odd and carry out the “‘integration’’ only 
from x4 = 3(a + b) toxy41 = 5. Then the correct result is obtained by replacing 
the elements Ly(Djx) and Ly —j41,n—e41(Dn-j41,n-241) both by their sum for 
j,& < M. In obtaining the numerical results reported in this paper this procedure 
was followed. The matrix a”,; was computed utilizing a numerical differentiating 
routine which uses polynomial approximations, but a procedure similar to that 
outlined for obtaining a°,; could have been followed. 

In both methods, the scheme used in step (2) for solving the equation 


Ly = \Dy 


begins with the computation of the characteristic values and characteristic vectors 
of the matrix D, to obtain the decomposition D = U7AU, where A is a diagonal 
matrix containing the characteristic values of D and U is a unitary matrix con- 
taining its characteristic vectors. Then 


[(A4U7)TL(A4U7) ](UAtY) = A(UAtY). 


The matrix shown in brackets is computed and then diagonalized to obtain 
the characteristic values 4 and the vectors UA'Y. The characteristic vectors 
Y = A*U7(UAtY) are then computed. Both diagonalizations are done by the 
Jacobi method (see Gregory [3 ]). However for Method I advantage is taken of the 
fact that D is already a diagonal matrix. 

Since, except for the factor rt in the definition of P*,,, the elements P*,,, Q*,:, 
R‘,, depend only on the sum r + ¢ and i, a total of (N + 1) (6m + 2) quadratures, 
each over an interval of length 4, are required. Computing the matrices L and D 
requires a total of m(N + 1)(N? + N + 2m) multiplications. The most trouble- 
some term in this expression is mN*; this is to be compared with 50N?*, which is 
approximately the number of multiplications required in the two diagonalizations. 
Thus the time required to compute the elements of L and D is small in comparison 
to the time required to solve the equation Ly = ADy. 
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10. Numerical Results and Their Discussion. In this section we shall report 
and discuss the results obtained by applying Methods I and II to the example 
given by equations (6.7) to (6.9). The results are given in tabular form. Each 
entry represents the value of 100(A,.* — Ax)/A, where , is the characteristic 
value as computed from equation (6.10) with = k and ,* is the kth approxi- 
mate characteristic value obtained by solving the equation 


Ly = Dy. 


Results on the comparison of the y; which satisfy the above equation and 
y(x;) computed from equation (6.11) are also available but not reported herein 
because of space limitations. 


TABLE 1. 100(Ax* — Ax)/Ax from Method I with N = 19 


k= 1 2 3 + 5 6 7 8 9 
ist Order 4.97 4.29 3.16 1.60 — 37 —2.75 —5.49 -—859  —11.99 
2nd Order 00 — 03* — 15 — 47* -1.17 -—2.46* -4.60 —-—7.89 —12.60 
2nd Order* .00 00* — 14 — 47* -1.17 -—245* -460 -—7.88  —12.59 


Not all the computed \*’s are reported in these tables. In case N = 19 we 
give nine \* and in case V = 9 we give five \*. The reason for this is that the true 
characteristic values \, for m > 10 when N = 19 and m > 5 when N = 9 corre- 
spond to functions which have a maximum within the first interval used. Hence 
any approximation of this function by the value of its ordinate at the end of the 
interval will not be a reasonable one and would be expected to lead to poor 
results, as was found to be the case. 

In Table 1 which contains the results obtained by the use of Method I with 
N = 19 still an additional selection principle was used in choosing the value of 
d* reported. The first nine values of \* which were associated with values of y;*, 


TABLE 2. 100(d,* — Xx)/Ax from Method II with N = 9 


m/k 1 2 3 4 5 
2 018 309 1.249 4.192 7.498 
3 001 052 645 2.224 7.299 
4 .000 015 022 323 1.945 
5 .000 .000 019 104 767 
6 .000 .000 .003 106 383 
7 .000 .000 .000 028 .723 
8 .000 .000 .000 .004 .203 


such that a continuous function passing through these ordinates would have the 
same number of zeros as the corresponding solution of the Sturm-Liouville 
equation, were used. That is, the known oscillation properties of solutions of 
Sturm-Liouville systems were used to distinguish between ‘‘true’’ approximations 
to \, and “‘spurious’’ ones when two close values of \,* were obtained. The entries 
marked with an asterisk in Table 1 were one of a pair of computed values of \* 
which were numerically close. 

The first line of Table 1 contains the results obtained by using the first order 
expressions for a;; and w; (equations (5.3) and (5.4) in equation (3.4)). The second 
line was obtained by using the second order expressions for a;; and w; (equations 
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(6.3) and (6.4)). The third line was obtained by modifying the even rows in the 
body of the matrix in equation (6.3) so that they are of the form 


2-16 _— (3 a 2-14) 0 (3 a 2-14) —2-16 
instead of 
o A 6.4 & 


The fact that the entries of the third row are uniformly smaller than those of 
the second shows that there is no advantage to be obtained in Method I by 
requiring the a;; and w; to be chosen so that the resulting algebraic equation is a 
second order approximation to the original differential equation. Indeed this 
example shows that by modifying the a;; and w; satisfying this condition, so that 
the condition is violated, better results can be obtained. 

Table 2 gives the results obtained by using Method II. The entries in a single 
row correspond to a fixed value of m labelling that row. As was pointed out earlier, 
the total computation time is not seriously changed as m is increased from 2 to 8. 
By using m = 8 one can obtain all the characteristic values one has a right to 
expect to get with an accuracy less than or equal to .203 per cent. 


TABLE 3. 100(Ax* — Ax)/Ax from Modified Method II (with m = 2) 


k= 1 2 3 4 5 6 7 8 9 
N=19 .00 —.01 — .05 — .13 — .32 —-67 -1.39 -3.09 —7.53 
N=9 —010 —.137 — 655 -—3.075 —18.943 
N=9 011 126 389 674 1.34 
N=9 .018 .308 1.249 4.192 7.498 


The first two rows in Table 3 contain results obtained by using a modification 
of Method II. In the first row of this table N = 19 and in the last three rows 
N = 9. The last row of Table 3 is the first row of Table 2. The modifications 
referred to are the following: The polynomials Y;(x) used to replace the functions 
y(x) in the integral (2.2) are such that 

Yo(x) = P2(x) 
Viga(%ig1) = Vi(xi4s) 
V" igs (%igs) = Ya’ (wi41). 


Further, for the first two rows of Table 3 the matrix D is made diagonal by 


evaluating 
6 
f ry*dx 


by means of Simpson’s rule instead of replacing y(x) by Y(x) defined above. 
This last approximation procedure means that the computation time required 
for solving the equation 


Ly = \Dy 


is markedly decreased. However, it is no longer true that the \,* are always 
greater than the d,. 
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On comparing the second and fourth rows of Table 3 we see that for the first 
four characteristic values the results obtained from the modified Method II are 
better than those given by the unmodified one for NV = 9. Other computations 
have shown that this result does not hold for m # 2 if the class of admissible 
function is defined so as to have continuous derivatives at the points x = x; and 
if the matrix D is evaluated as in Method II. 

In the third row of Table 3 the matrix D is evaluated as in Method II where 
the admissible functions are defined as for the first two rows. 

The numerical results, although obtained only for a particular simple example, 
would seem to justify the use of the methods discussed above for the solution of 
more complicated Sturm-Liouville problems. They should be supplemented by a 
test for the accuracy of the results obtained for the characteristic values and 
characteristic functions applicable to problems for which the solution is not 
known. Various such tests are under investigation and will be reported on 
subsequently. 
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On the Numerical Evaluation of the Stokes’ 
Stream Function 


1. Introduction. In the study of axially symmetric problems in fluid dynamics, 
the Stokes’ stream function, that is, the function which satisfies 


1 
(1.1) tan Ste tthe eM y ~ 0, 


is of considerabic interest. This function is constant on the streamlines. 

The problem to be considered is a Dirichlet type problem. Let G be a closed, 
bounded, simply connected plane region whose interior is denoted by R and 
whose boundary curve is denoted by S. Let G not contain any point where y = 0. 
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Let¥g (x, y) be defined and continuous on S. The problem then is to produce a 
function «(x, y) such that 


a) u(x, y) = g(x, y), on S, and 
b) u(x, y) satisfies (1.1) in R. 


Under general conditions, there exists a unique solution (Bernstein [1] and 
Courant and Hilbert [3°]), and only such cases will be considered. However, the 
analytical determination of u(x, y) is quite another story from that of its existence 
and usually offers what are at present insurmountable problems. The approach 
here, then, will be from a numerical analysis point of view. 


Po ee 
kf i... 


\ 5 6 7 8 9 


By 11 12 13 | 
15 16 17 a, 


{My 21 22 wo, 


DracraM 1 
























































Let h be a fixed positive constant called the mesh size. Let (xo, yo) be an 
arbitrary, but fixed point of G and denote by G, the set of all points of form 
(xo + mh, yo + nh), contained in G, where m and n are integers. Two points (x1, 1) 
and (x2, y2) of G are called adjacent if and only if (x2 — x1)? + (ye — y:)? = F#*. 
R,, the interior of G,, is the set of all points of G, whose four adjacent points 
belong to Gy. S,, the boundary of G;, called the lattice boundary, is defined by 
S; = G, — R,. It is also assumed that any pair of points of R, can be joined by a 
connected polygonal arc consisting of straight line segments which join adjacent 
points of R,. 

The technique proposed will involve the application of a difference equation 
to yield a system of linear equations whose solution is an approximation to the 
solution of the formulated Dirichlet problem at the points of Gy. 

The question of solving the linear system will not be considered. Systems with 
unique solutions may be approached by means of Cramer’s rule, Gauss’ Elimina- 
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tion Procedure, matrix inversion, relaxation, iteration, gradient methods, and 
other techniques (Geiringer [4], Greenspan [6], Milne [7], Shortley et al [8], 
NBS [9], Young [13, 14]). 

The author is indebted to D. M. Young for valuable criticisms and suggestions. 

2. General Method. Suppose G, consists of points. Number these points 
in a one-to-one fashion with the integers 1, 2, 3, ---, . Denote the coordinates 
of the point numbered k by (x:, y,) and the unknown stream function at (xx, ys) 
by (xz, yx) = ux, for k = 1,2, ---, n. 

Let (x;, yi) be an arbitrary point of S,, the lattice boundary. Approximate u; 
by g(x’, y’), where (x’, y’) is the nearest point of S to (x;, y,). If (x’, y’) is not 
unique, choose any one of the set of nearest points and use it. The problem of 








n y 
2 _ 
1 
—}- of + >x 
3 (x,y) 
on ae 4 
DIAGRAM 2 


finding numerical approximations to u(x, y) on the lattice boundary is, though 
crudely done, adequate for present purposes. For diagram 1 this means that at 
the points 1, 2, 3, 4, 5, 7, 9, 10, 14, 15, 19, 20, 21, 22, 23 the values of the stream 
function have been approximated by values at nearby points of S. 

We then require that at each point (x;, y;) of R,, the function u satisfy 


(2.1) 4u(x;, yi) — u(xi + h, yi) — (1 _ +) u(xi,yi th) ~ 
— u(x; — h, yi) — (1 _ F) wie, yi — h) = 0. 
24: 


It is convenient, in practice, to use the subscript notation. For example, for 
diagram 1, if 7 = 12, (2.1) becomes 


h h 
Auy. — U13 — (1 ait rae - a, (1 +t) un = 0. 
2y12 2y12 
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Note that w12, 413, 411, 417 are unknowns while uw; is a constant determined by the 
method for points of S,. Application of (2.1) to each point of R, thereby results 
in a system of linear equations which, when solved, yields the remaining numerical 
approximations. 

Equation (2.1) is only one of a variety of possible difference equations which 
can be developed, as described in the next section. 

3. Development of a General Difference Equation. We shall seek a solution 
u(x, y) which is of class C* and use the notation 0*+*u/dx‘dy’ = u;;,i + 7 < 4, 
and Ai = Us, ;/t!7!. 

Let (x, y) be an arbitrary point of R,. The aim will be to develop a “‘5-point 
formula,” so consider here only (x,y) and its four adjacent points. Let (x,y), 
(x +h, y), (x,y +h), (xn — h,y), and (x,y — h) be denoted, respectively, by 
0, 1, 2, 3, 4, as in diagram 2. Hence (xo, yo) = (x,y), (x1, 91) = (x +h, y), 
(x2, y2) = (x,y +h), (xs, ys) = (x — hy), (x4, 94) = (x, » — A). Also let O(h") 
represent any and all functions of order at least m in h. 

In order to deduce a difference equation which approximates (1.1), let 


4 
(3.1) L(u) = ¥ any. 
0 


By use of the finite Taylor series expansions about (xo, yo) for the various %; 
and by the definition of A;;, (3.1) becomes 


(3.2) L(u) = Aoo(ao + a1 + a2 + a3 + a) + hLA 10(a1 _ a3) + Aoi(a2 - as) | 
+ W[A2(a1 + as) + Ao2(az + as) ] 

+ h®CA g0(a1 sais a3) + Ao3(a2 ™ as) | + 0. (h*) 

Equation (1.1) may be rewritten as yuoo — uo,1 + yuo,2 = 0. Since u is 


assumed to be of class C‘, further differentiation of (1.1) yields 
Yus,o — U1,1 + YUr,2 = 0 
Uo,o + Yuo,1 + yuo,s = 0. 

Rewriting these in terms of the A;;’s and solving, one finds 


Ao Au Aw Ao Ag 
33) An =—"*-—Aai An =" =="; An = — 34a 4+—" —-—. 
( ) 20 2y 02 30 6y 3 21 03 y 2 


Substitution of (3.3) into (3.2) yields 
(3.4) L(u) = Aoo(ao + a1 + a2 + a3 + a4) 
h 
+h [40 (= + = + az — a) + Axo(a: — a) | 
2y ay 
+ h®[Aos(a2 + a4 — a1 — a) J 
+ h®[Aso(a1 — as) + Aos(a2 — a4) ] + O(h'). 
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Let 
ao + a1 + ae + 3 + oy = 1 
h 
(3.5) — oR 4 aa" @& 
2y say 
ai as = €3 
a2 + a4 — a1 — a3 = & 
a2— a = €5, 
where 


a. = Off), e = O(), « = O(h*), « = Oh), & = Ofh). 


The solution to system (3.5) is 


-97_ 9,8 Wa ei 9 
aT KP Te mmapigeg ss 
~7 98 & | _ Wei w_ 

(3.6) iat ek et Wile tat it foe: 


a =a — a t4 — 4%. 


Finally, setting (3.1) equal to zero and substituting (3.6) in (3.1) yield the general 
difference equation (3.7) which approximates the differential equation (1.1) 


(3.7) (« - uo ae - 2) 4 +(2 ~ 242) a, 


hook .". 
a 4 ey ey a 4 
+(2+H+5 2) u+(% h ) a 
ey a «a _ 
+(2+8 2 a . 


The process, then, was to find a; such that: L(u) + Ofh*) = 0, where x satisfies 
(1.1), and to then use Z(u) = 0 for the general approximation. 

Equation (2.1), which will be used throughout, is deduced from (3.7) by letting 
€1 = €2 = €3 = eg = 0, es = h, by dividing through by y, and by noting u» = u(x, y), 
uy = u(x +h, y), ue = u(x, y +h), us = u(x — hh, y), Ug = U(x, y — A). 

Other difference equations may be constructed by selecting the e,’s differently. 
One of these of particular interest is one which would yield a linear system which 
has a symmetric coefficient matrix. This may be done in the following manner. 
Number the points as in diagram 1 so that the numbers increase from left to 
right and also increase if one reads down any column of points. Apply difference 
equation 

2 2 2 Uy U3 
(3.8) w(+24572 74575) . 


ie on = 2 h 
m( U4 2y or 0, | y | # 
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to the points of R,. Start with the first, or highest, row of points; traverse the 
row from left to right; proceed to the next row; traverse the row from left to 
right; proceed to the next row, etc.; continue the process until all the points of 
R, have been used. Equation (3.8) may be constructed by a numerical procedure 
for certain self-adjoint elliptic differential equations (Young [15]), or directly 
from (3.7) by letting 











0 — + = 0 
= . => 9 = . 
. eo” ¥ 7 Qy+hQy—h)’ © 
—8y re 2 +4h 
€ = watt € = . 
* Qy +h)(Qy—h) " y’? * (2y + h)(2y — h) 


4. A Theorem on an Error Bound. In this section, let u(x, y) be the solution 
of the Dirichlet problem being considered and U(x, y) the solution of the nu- 
merical method described in section 2. This means that 


h h 
(4.1) -4U+0itut+(1-2)u+(1++)u-0 
2y 2y 


for any arrangement of points of G; like that of diagram 2. 
Lemma 1. The solution of the system of linear equations which results by applica- 


tion of the process of section 2 with y ¥ 0 and 


< 1 ts unique. 








Proof. It is sufficient to show that the determinant of the system of linear 
equations is not zero and this is done by demonstrating that the only solution of 
the homogeneous system which results by considering g(x,y) = 0 on S is the 
zero vector. Suppose then there exists a nontrivial solution for the homogeneous 
system. For some point of R,, U # 0. Suppose U > 0. Let the largest value M 
occur at (xo, yo). Then 


h h 

(4.2) 4U, = Ui + u+(1-*)n+(14+4)u 
2y 2y 

and 

(4.3) M=U.> Ui, i=1,2,3,4. 


Case 1. U2 = U,. Then (4.1) becomes 
(4.4) 4U,) = Ui + U2+ Us + Ui. 


But (4.3) and (4.4) imply Up = Ui; = U2 = U; = Uy, = M. Hence one can 
apply the same argument at (x:, yi) and continue in a finite number of steps to 
show that U at a boundary point is equal to M. This is a contradiction since 
g(x, y) = 0 on Sand the method described for points of S, implies U at all lattice 
boundary points is zero. 
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Case 2. Uz > Uy. Let Us = U2 — kh, k > 0. Then 


h 
4U, = U, + 2U2+ v,-e(1+2), 
and 


h 
4Uy < U, + 2U2+ U3, since ll lia 


But this is impossible since U; < Uo, U2 < Uo, Uz < Us. 

Case 3. Us > U2. The discussion is analogous to that for Case 2 and leads to 
a contradiction. 

Hence, in all three cases, the assumption that U > 0 yields a contradiction. 
Contradictions are similarly reached if U < 0 is assumed. Hence the only solution 
is the trivial solution. This completes the proof. 

Using the same ideas and techniques first set forth by Gershgorin [5], we 
apply these by assuming in all that follows that ¢ = h/2y, |o|< 1 and by defining 


(4.5) ALu(x, y)] = hk? — 4u(x, y) + u(x +h, y) + u(x — h, y) 


h h 
+(1 —F) uy +0 + (14%) aces - mf. 


LEMMA 2. If A[v] < 0 on R, and v > 0 on Sj, then v > 0 on Ry. 
Proof. Suppose v < 0 for some point of R,. Then there exists a point (xo, yo) 
of R, where v is a minimum. Hence ° 


(4.6) v(x0, Yo) < 0, and, v(xo, yo) < v(x, y), for all (x, y) of Ry. 


Now, A[v] < 0, by assumption, so that, equivalently 


1 h h 
(4.7) w>ifn+(1-2)mtu+(1+4) al. 


But (4.6) and (4.7) imply, in a manner similar to that used in Lemma 1, 
that v9 = v1; = v2 = vs = v. In a finite number of steps this leads to the result 
that v < 0 for some point of S;, which is a contradiction. This completes the proof. 

Lemma 3. If —|A[v1]| > Ave] im R, and |v1| < v2 om Sh, thén |v1| < v2 on Ry. 

Proof. v2 —-%>0 on S, and Ave — v1] = ALve] — AL]. Since Alve] 
+|A(v]| < 0, then A[ve] — Afi] < Ale] +|AL1]| < 0. Hence ALv2 — v1] < 0. 
By Lemma 2, v2 — v1 > 0 on R,. 

Also, ve + v1 > ve —|01|> 0 on S,. Then, ALve + v1] = ALve] + Af] < 0. 
Hence v2 + 7; > 0 on Ry. But ve + v7; > 0 on R, and ve — v1; > 0 on R, imply 
v2 >|v1| on Ry. This completes the proof. 

Lemma 4. If |A[v]| <A on R, and if |v| <B on Sy, and r is the radius of a 


A 
properly selected circle, as described below, which contains G, then |\v| < a r+ B. 
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Proof. Let 


w(x, 9) = [4 {1 - ea ait o— wl + al, 





where (x — xo)? + (y — yo)? = r° is the equation of any circle containing G for 
which yo/y > 1, for all yin G. At least one such circle exists since G is simply 
connected, bounded, closed and y ¥ 0, which imply that y is of constant sign. 


A 
Direct calculation yields A[w] = — 2 (1+ 0/y) < —A. Also, w > Bon S;. 


Now, since |A[v]|< A on R,, by assumption, and it has been shown that 
A[w] < —A, or equivalently, — A[w] > A, it follows that — A[w] > — A > Av] 
on R,. Since |v| < Bon S;,andw > Bon S;,w >|v|. Hence, — |ALv]| >ALw] 


Ar 
on R, and w >|v| on Ss. By Lemma 3, w >|»| on R,, or |v] < w < = +B. 


This completes the proof. 
THEOREM. If u(x, y) is of class C® in a closed region G, y # 0 in G, u denotes 
the solution of the stream equation, U denotes the solution of the linear system which 














h 
results from application of (2.1) in the method of section 2, |o| = | < 1, then 
¥y 
W?M; . h®?M, | 
i 2hM;, 
bd ui so Si + 32 7 ee 
where 
M; = Max [Max |1,0|, Max|wo1|], Ms = Max|uos|, 
(z,y)G@ (z,y) eG (z,y) eG 
M,= Max [ Max | wo,4|, Max | #4 0] J, M; = Max | u,5|, 


(z,y) eG (z,y) G (z,y)G 
7 = GLB\y| in G, r = radius of any circle of type described in Lemma 4. 
1 —_— : 
Proof. Let R = \[u] — (uw. - = + uy } . Substitution of the finite Taylor 


series expansions for “1, 2, 3, % in AL“ ] and use of the Mean Value Theorem 
yields 


h? 
R= rT [4 0(01, ¥) + to, a(x, p2) + ta, 0(ps, ¥) + to, 4(%, pa) ] 
! a 
= Bly Uae S¢ = 2. 75 [uto, s(x, ps) J, 
where p2 < ps < ps. Hence 


Pt ew 4 Mist 
3! 319 | 2-449 





IR|< 


1 
Since cn Be Uy, = 0, 


1 2M, 382M, kM, 
ani == <« = < — - 
(4.8) or (wa * uy + t)| |A[u]| =|R| < —— 349 4 31 + oa 
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Also, for any point of S,, U was taken as the value of g(x’, y’) at the nearest 
point (x’, y’) on the boundary S, and g(x’, y’) = u(x’, y’) on S. Thereby, for any 
point (x, y) of Sh, 
(4.9) |U(x,y) — u(x, y)|=|e(x’, »’) — u(x, »)| =|u(x’, »’) — u(x, y)|, 
where (x’, y’) is a point of S and (x — x’)? + (y — y’)? < 2h*. Therefore 


(4.10) |U(x,y) — u(x, y)|=|u(e’, y’) — u(x, y’) + u(x, 9’) — u(x, »)| 
S | he (p0, 9") (@ — x")| + | my (x, 02) (y — 9')| 
<MiC\y’ — y|+|x’ — x|] < 20M,. 


It must also be noted that ATU] = 0, by (4.1). Hence 
(4.11) |A[w — U]| =|ALu] — ACU]| =|ALu]]. 


Applying Lemma 4 to (4.8), (4.10), and (4.11), one finds that on R, 











2 2M 2 4 
“\t 3 hW®M, AM, + 2h. 


-al<= 
Gm lU- 439) 39 + ar tray 


This completes the proof. 


CorROLLARY. Under the assumptions of the theorem just proved, U > u,ash > 0. 

Under general conditions set forth by Collatz [2], the term 24M, in the error 
bound for |U — u| may be replaced by a term of type O(h*) by replacing the 
unesthetic method for deciding U on S, as follows. Let (x1, y:) be a point of S, 
(x2, y2) of Ry, and (xo, yo) of S,. Let PQ (diagram 3) be parallel to a coordinate 
axis and let 6 be the absolute value of the distance between (xo, yo) and (x1, ¥1). 
Choose Up by 


_ Ag (es, 91) + 602 
h+é 





Uo 


5. Use of Diagonal Points. After a numerical solution has been found on a 
given set G,, solutions at other points may be approximated rapidly as follows. 
Let 5, 6, 7, 8 represent four points of G, which are vertices of a square. (See 
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diagram 4.) Let 0 represent the diagonal point of the square. Use 


(2 — o — o?)(Us + Us) + (2+ 0 — o°)(Ur + Us) J 
(5.1) Us= 4(2 — o?) —_" 


The unexciting calculations will here be replaced by a general outline of how to 
establish (5.1). Let L(u) = aotuo + asus + acttg + arttz + agtg. Expand uss, 
uz, us in Taylor series about (xo, yo) and substitute these expressions into L(u). 
Simplify by replacing u,;, ; by 7!7!A.;. Apply identities (4.3) and combine as in (4.4). 
Set the coefficients equal to ¢;’s. This system of equations is 


ao + as + as + a7 + ag = €1, — a — ar + a3 = €2, a5 + os — a7 — Og = €5, 


a5 
cs s @€2 
u(1+5-$)+a(1+5-$) 


o o 
tar(-14£4+5) +a(-1454+5) =0, 


w(+4)-a(049) 40-2) -a(-2) = 


where e; = O(h*), 2 = O(h*), es = O(h®), eg = O(h*), «5 = O(h). A general formula 
can be established, but by letting «1 = «2 = ¢: = e« = 0, es = h, solving for the 
a,’s and using the method of section 3, one deduces (5.1). 

6. Concluding Remarks. The method and ideas developed can be used in an 
analogous fashion to effect the numerical solution of problems involving any one 
of the class of equations 





(6.1) Ure + My + ty = 0, 
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where & is a fixed constant. Such equations are confronted frequently in physics 
(Shortley et al [8]). 

Finally, two problems remain which are worthy of note. Many physical situa- 
tions allow y to be zero on the boundary S and a solution which satisfies (6.1) 
in R is desired. An effective technique in the sense that U — u as h — 0 is needed. 
Also, a reasonable method for approximating the error in the numerical solution 
is wanting. The error bound (4.12) involves unknown quantities and the only 
work which has been done toward establishing error bounds which can be calcu- 
lated has been for harmonic functions (Walsh and Young [11, 121). 


DoNALD GREENSPAN 
Hughes Aircraft Company 
Culver City, California 
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Tables of Values of 16 Integrals of Algebraic- 
Hyperbolic Type 


This paper gives tables of values of the following 16 integrals of algebraic- 
hyperbolic type. 


' bi... =f xkdx (k > 1) 
(1) T* k!Jo sinh2x+2x’ (k > 3) 
. I, _ 2 2c xte-2=dx (k > 1) 
(2) Ih*  kiJo sinh 2x + 2x’ (k > 3) 
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(3) i _ 2 (* x* tanh xdx (k > 0) 
Il,* k!Jo sinh 2x + 2x’ (k > 2) 
(4) IV, _ 2* (* x* coth xdx (k > 2) 
IV.* kiJo sinh2x+2x’ (k > 4) 
5) Vi wee px x* sinh xdx (k > 0) 
Vit kiJo sinh2x+2x’ (k > 2) 
6) Vin _ © a. x* cosh xdx (k > 1) 
VIi.* kiJo sinh2x+2x’ (k > 3) 
(7) VII, _ py x*tanhxsinhxdx (k > 0) 
VI* kiJo sinh2x+2x ’ (k>1) 
8) VII, _ 1 @“x* cothxcoshxdx (k > 2) 
VIITL,* k!Jo — sinh 2x + 2x (k > 4) 
The integral (7) also converges for k = —1. It can be shown that 
° tanh x sinh xdx * 1— IIe 
=In3 — ———— = 0. 7. 
J, x(sinh 2x + 2x) ' 2 (2m + 1)2?" a 


With the factors as shown, the integrals tend to unity as the integer k tends 
to infinity, except JJ, and JJ,* both of which tend asymptotically to 2-@. 
The first four integrals, which are called Howland’s integrals, arise from the 
problem of a symmetrically perforated strip or a notched strip (Howland [1], 


TABLE 1 
k Ik Ii, Ili; IV: Vi Vik Vil; VII, 
0 a= — 0.29662 0 — 0.52685 6 —_— 0.36980 1 -— 
1 0.76857 5 0.22012 0 0.47443 0 — 0.73884 4 0.91747 7 0.65893 5 = 
2 0.76784 7 0.08792 7 0.63084 1 1.09596 7 0.86986 6 0.92259 0 0.83403 0 1.02797 8 
3 0.82771 0 0.04334 8 0.75383 6 0.94313 8 0.94005 4 0.95831 1 0.92528 4 0.98312 4 
4 0.88350 7 0.02258 3 0.84283 9 0.93632 3 0.97401 0 0.98049 8 0.96827 5 0.98810 2 
5 0.92547 6 0.01192 3 0.90323 7 0.95164 4 0.98924 9 0.99154 5 0.98711 6 0.99405 4 
6 0.95419 2 0.00628 8 0.94219 2 0.96751 5 0.99571 1 0.99651 5 0.99494 2 0.99736 1 
7 0.97269 9 0.00329 5 0.96631 5 9.97953 8 0.99833 7 0.99861 5 0.99806 6 0.99890 2 
8 0.98412 4 0.00171 3 0.98077 1 0.98763 5 0.99936 9 0.99946 4 0.99927 5 0.99956 1 
9 0.99094 9 0.00088 4 0.98920 7 0.99274 6 0.99976 5 0.99979 7 0.99973 2 0.99983 0 
10 0.99492 2 0.00045 3 0.99402 5 9.99583 9 0.99991 3 0.99992 4 0.99990 2 0.99993 5 
11 0.99718 9 0.00023 1 0.99673 0 0.99765 4 0.99996 8 0.99997 2 0.99996 5 0.99997 6 
12 0.99846 0 0.00011 7 0.99822 7 0.99869 5 0.99998 9 0.99999 0 0.99998 7 0.99999 1 
13 0.99916 4 0.00005 9 0.99904 6 0.99928 3 0.99999 6 0.99999 6 0.99999 5 0.99999 7 
14 0.99954 9 0.00003 0 0.99949 0 0.99960 9 0.99999 9 0.99999 9 0.99999 8 0.99999 9 
15 0.99975 9 0.00001 5 0.99972 9 0.99978 9 0.99999 9 1.00000 0 0.99999 9 1.00000 0 
16 0.99987 1 0.00000 8 0.99985 6 0.99988 6 1.00000 0 1.00000 0 
17 0.99993 2 0.00000 4 0.99992 4 0.99993 9 
18 0.99996 4 0.00000 2 0.99996 0 0.99996 8 
19 0.99998 1 0.00000 1 0.99997 9 0.99998 3 
20 0.99999 0 0.00000 0 0.99998 9 0.99999 1 
21 0.99999 5 0.99999 4 0.99999 5 
22 0.99999 7 0.99999 7 0.99999 7 
23 0.99999 9 0.99999 8 0.99999 9 
24 0.99999 9 0.99999 9 0.99999 9 
25 1.00000 0 1.00000 0 1.00000 0 
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TABLE 2 

k is II,* ITI;* IV;* V;* Vi* VII,* VIII;* 
1 _ _ _ _ _ _ 1.51115 5 _ 

2 _ _ 2.13561 8 —_ 1.40879 6 —_ 1.15546 4 _— 

3 2.03871 1 0.46071 4 1.41506 3 — 1.10522 6 1.20400 4 1.05738 7 — 

4 1.35329 4 0.09931 6 1.19555 3 1.70756 9 1.03438 8 1.05121 3 1.02214 3 1.08171 9 
5 1.15686 4 0.03241 3 1.10049 3 1.24000 1 1.01212 3 1.01624 8 1.00860 2 1.02146 1 
6 1.07673 0 0.01261 7 1.05358 0 1.10538 5 1.00439 8 1.00556 4 1.00332 5 1.00686 4 
7 1.03925 1 0.00539 1 1.02903 4 1.05082 2 1.00161 1 1.00196 5 1.00127 3 1.00233 9 
8 1.02053 8 0.00243 3 1.01583 4 1.02560 8 1.00059 1 1.00070 3 1.00048 3 1.00081 7 
9 1.01087 0 0.00113 6 1.00864 7 1.01319 9 1.00021 6 1.00025 2 1.00018 1 1.00028 9 
10 1.00578 5 0.00054 2 1.00471 5 1.00688 6 1.00007 9 1.00009 1 1.00006 7 1.00010 2 
11 1.00308 5 0.00026 3 1.00256 4 1.00361 5 1.00002 9 1.00003 2 1.00002 5 1.00003 6 
12 1.00164 5 0.00012 8 1.00139 0 1.00190 3 1.00001 0 1.00001 2 1.00000 9 1.00001 3 
13 1.00087 6 0.00006 3 1.00075 0 1.00100 3 1.00000 4 1.00000 4 1.00000 3 1.00000 5 
14 1.00046 6 0.00003 1 1.00040 3 1.00052 9 1.00000 1 1.00000 1 1.00000 1 1.00000 2 
15 1.00024 7 0.00001 6 1.00021 6 1.00027 8 1.00000 0 1.00000 1 1.00000 0 1.00000 1 
16 1.00013 1 0.00000 8 1.00011 5 1.00014 6 1.00000 0 1.00000 0 
17 1.00006 9 0.00000 4 1.00006 1 1.00007 7 

18 1.00003 6 0.00000 2 1.00003 3 1.00004 0 

19 1.00001 9 0.00000 1 1.00001 7 1.00002 1 
20 1.00001 0 0.00000 0 1.00000 9 1.00001 1 
21 1.00000 5 1.00000 5 1.00000 6 
22 1.00000 3 1.00000 3 1.00000 3 

23 1.00000 1 1.00000 1 1.00000 2 
24 1.00000 1 1.00000 1 1.00000 1 
25 1.00000. 0 1.00000 0 1.00000 0 


Howland and Stevenson [2], Ling [3, 4]). The remaining twelve integrals, save 
VII, and VIII,*, arise from the problem of an unsymmetrically perforated strip 
(Ling [5 ]). 

The four Howland integrals J,, I,*, IJ;, I[;,* were tabulated by Howland [1] 
and Howland and Stevenson [2] to 5D. In a recent paper by Nelson and the 
present writer [6] these integrals were recalculated and tabulated to 6D. The 
values of J, when k is odd were also tabulated by Nelson [7] to 9D. (However 
our J; is designated by Nelson as a;. Nelson has also calculated the remaining 
Howland integrals to 9D, but this work is unpublished.) For convenience of 
reference, values of the four Howland integrals are reproduced in Tables 1 and 2. 
It will be shown that the remaining twelve integrals can be evaluated in terms of 
the four Howland integrals. 

The four integrals I7J,, JJI,*, IV, IVi* may be evaluated by splitting the 
integrands and then integrating from zero to infinity, 








x*tanhx _ xt i x*-1(1 — e-*) ah—le-2 
(0) sinh 2x + 2x sinh 2x + 2x 2(sinh 2x + 2x)  2coshx’ 
x*cothx _ x x*-1(1 + e-**) xi-ig- 








sinh 2x + 2x sinh 2x + 2x 2(sinh 2x + 2x) ~ 2 sinh x 
The integrals become 


TTT, = In + Uae — Una — Sx)/k, 


(10) ITi,* = I,* — (I*,-1 7 ITIT* 1 = Sx)/k, 
IVe = Te — e-s + ITn-1 — Sx)/k, 
IV»* = Te* + (I*e-1 + TT — Si)/k, 
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where 
Ren Qk-1 ft x*-le-zdx . 1 1 4 1 1 + 
an)  @= tt coshx ae" 3h 4h : 
aid Qe-1 foe 214545454 
*(k—-1)!Jo sinhx ae * ge * ge 


The values of s, and S, were tabulated by Glaisher [8 ]. 
Next, consider the four integrals Vi, Vi*, Vik, VJ.*. By expanding sinh x 
into power series of x, the integral V; develops into the series 





2n + k + ‘) Tonpe+t 


Jrntk+i ‘ 


(12) Vi = E ( 


n=O B 


In order to improve the convergence, Kummer transformation (Knopp [9]) may 
be used. The series then becomes 


1 ms 2n+k+1 1 — Tengesi 
(13) Me= 1-55 -E( k =. 


Bet 


Similarly, the integral V,* develops into the series 


2n a k _ 1) Fomtecs — 1 


k 22ntk+i 


(14) Vit =1- get E( 


By expanding cosh x into power series of x, the following series are obtained 
in the same way, 


2n + k Tonite 
VE = 1455-2 ("F") ie. 
2n a k I* one —1 
k Qentk 


(15) 


VI.* lt+gntt 


gen 


The last four integrals may be evaluated by splitting the integrands into the 
following and then integrating from zero to infinity, 





x* tanh x sinh x 7 x* cosh x x*—! sinh x x*-1 
(16) sinh 2x + 2x  sinh2x+2x  sinh2x+2x 2coshx’ 
x* coth x cosh x x* sinh x x*—! cosh x el 


= B= ‘ 
sinh 2x + 2x sinh 2x + 2x sinh 2x + 2x = 2 sinh x 





The integrals become 
VII, = VIk + (Vi-1 — ue)/k, 


ai VIT.* = VIp* — (V*e-1 — m)/k, 
VITI, = Vi _ (VIp-1 _ U;)/k, 
VIIL,* = Vi* + (VI*R-1 — U:)/k, 
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where 
nit 1 ee ele oN 
ions * " 2(k —1)!J0 coshx ge" 5 7k ' 
1 f= 1 1 1 
= =1 - os — ee, 
Us 2(k —1)!J0 sinh x ta tyt at 


The values of U; and “ were also tabulated by Glaisher [8, 10]. 
The preceding method of evaluation is satisfactory as far as it goes, except that it 
does not cover the initial integrals J7Io, III,, III.*, III;*, VIIo, VII,*, and VII". 
The integrals J7I) and III, may be evaluated by expanding tanh x into the 
following series and then integrating from zero to infinity, 
1 = (3s)* 
(19) ae sinh 2x (-, (2m)! 


This leads to 





IIIg = & (Uen — Ton—1)/2n, 
n=1 


(20) 
ITT, = Y (Uen41 — Ion). 
n=1 
An alternative expression for JJ, may be obtained as follows. By using the 
expansion 





La) 2n 
(21) sinh x = tanh x - <r : 
the integral V; develops into the series 
1 Oo f[antk\1 — Fon 
(22) Vi = 1+gn- 2 ( k ja efits ? 
which gives when k = 1, 
7 2 2n+1 
(23) III, = ser 24-Vj)+> pe (1 — TDI on41). 


n=1 


The integrals JJJ;* and JII;* may be evaluated from the following series 
which is obtained in the same way, 
1 © ( Wn+k\ITT* one — 1 
(24) We mt gat E( Met eat ee 


When k = 2 and 3 respectively, this series gives 


23 iy ie 2 *ont2 — 1 
mre = 3 4 ave — 1) - E( n+ jf) 


ont 2 Q2n ? 
“) 73 © {2n+3\III* 1 
SES on ae + OCF — 0) at 
Int = + 8(V* — 1) - 5 ( : ) = 


n=1 


By expanding sinh x into power series of x, the integral VII) develops into 
the series 
2 - 1— III 
(26) VIIy = ;> pa et 


2n+1 
n=0 2" 





(3 


le 





ne 
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Similarly, the integral VJI,,* develops into the series 
2n + k + 1 TIT* on 4441 —1 
Qintkti ‘ 


VIkt =1— a5 +E k 


(7 zr 
When & = 1 and 2 Tr it gives 
> 7 
VIL,* = = = + x pa (EI* — 1), 
(28) 2n(2n + 1) 
© 2n(2n 
VII.* = —  o » rr) (IIT* n41 — 1). 


= 
Values of the twelve integrals thus computed and rounded off to six decimals 


iusng Nelson’s tables [7 ] when necessary, are given in Tables 1 and 2. The follow- 
ngi formulas are useful for checking purposes: 
214 — Nh) + SY CT T*¥ xy — [1 Ty1) = 2 

k=1 


eT 
(i — III) + & 2 (III*». — IIIx) = 1, 


k=2 


(1 — IVs) + E 5 (Vn — Vx) =1 
k=2 
@k(2k + 1 
(1 = 1) + EF (vgs — TV ngs) =1 


k=2 


o k(VI — Vi) = 3h, 


k=1 
ya (2k + 1) (VI ox41 = VITx41) = 4], 
(29) 
es DL k(VITIx — Vu) = ti, 
k=l 
= k(2k — 1)(VIu — Vilu) = Hs, 


k=1 


© k(2k +1 
ee ere, ~ Fars © dele 





a3 
5 ha = ve tr ar en 
= 
to x ASS (ve — vin) = HH, 
B(k = 1b — 208 = 3) ype — yt) = al Ve 





z 24 
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Tables for the Rapid and Accurate Numerical 
Evaluation of Certain Infinite Integrals 
Involving Bessel Functions 


Introduction. In a recent paper [1 ] the author has formulated a method based 
on Euler’s transformation of slowly convergent alternating series for the numerical 
evaluation of integrals of the form /.* f(x)dx, where a is a constant and where 
f(x) oscillates about zero in such a way that the integral over each half-cycle is 
smaller in absolute magnitude than (and opposite in sign to) that over the 
preceding half-cycle. The author has had occasion to make much use of this 
method for the evaluation of integrals of the type 


(1) [> rear, 
and ; 
(2) f * Ta(x)h(x)dex, 


where g(x), (x) are well-behaved continuous functions which tend to a finite 
constant value or zero as x tends to infinity ; here Jo(x), J:1(x) are Bessel Functions 
(of the first kind) of orders zero and one, respectively. The present paper gives 
tables useful in the evaluation of (1), (2). 

The author is grateful to Yigal Accad who performed most of the numerical 
calculations for Tables 1 and 2. 
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Description of the Method. The method employed [1] in evaluating integrals 
of the type (1) involved performing the integration over each of the first twenty 
half-cycles, i.e., evaluating the integrals 


(3) 2 Jo(x)e(x)dx i= 1,2, ---,20 


where xo is zero and x; is the ith zero of Jo(x). The first twenty terms of a slowly 
convergent alternating series for (1) were thus obtained. Euler’s transformation 
(Bromwich [2], p. 62) was then applied to this series in order to obtain a rapidly 
convergent series for the numerical value of (1). Integrals of the type (2) were 
treated similarly. 

In order to obtain high accuracy in the separate integrations (3) the Gauss 
quadrature formula (NBS AMS 37 [3]) was used for sixteen points of sub- 
division of each interval. Tables of abscissae and coefficients are given in [3] 
for nm = 2(1)16, 15D, and in Davis and Rabinowitz [7] for » = 2, 4, 8, 16, 20, 
24, 32, 40, 48, 20D. Values are also available in [7] for » = 64, 80 and 96. Fora 
given value of , Gauss’ formula provides an approximation equivalent to re- 
placing the integrand by a polynomial of degree 2m — 1. 

Description and Use of the Tables. Table 1 gives values of x and the corre- 
sponding values of Jo(x). The values of x were obtained by computation from 
values given in [3] for the interval (—1, 1) by means of the formula 


12—? , at? 
Te, + 2 





x =z 


where ~, g are the lower and upper limits of integration, and in our case are zeros 
of Jo(x) (the first value of p is zero), and the x; (given in [3 ]) refer to the interval 
(—1, 1). The zeros of Jo(x) are also given. 

The values of Jo(x) were obtained by interpolation from the Harvard Uni- 
versity tables [4]. In view of the high accuracy needed interpolation was effected 
by means of the first four terms of Taylor’s theorem 


he hi 
Jo(x + h) = Jo(x) + Jo (x) + 5 Jo" (x) + 5 Jo") 


- Ia(e) ) 


= Jo(x) = hJ1(x) +e yw (2 





+ oe (J +22 — 22@)), 
% oe 
these terms being sufficient for 10 place accuracy with the maximum value of h 
required. The tables are believed to be accurate to within a few units in the tenth 
place of decimals. 
Table 2 was similarly prepared for J:(x). Here the interpolation for J;(x) 
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was effected from [4] by means of the Taylor’s theorem expansion 


x66) ) + ye (22) be Jo(x) 
x x x 





Ji(x + h) = I(x) + &( Jute) - ~ r2)) 


3Jo(x) 65i(x) . 2Si(x) 
2 )©=6 + x - Jz). 





+1 ( 


In Tables 1 and 2 the zeros of Jo(x), Ji(x) were taken from tables in 
BAASMTC [5]. 
Table 3 gives the appropriate integration coefficients (these are quoted from 


TABLE 3. Gauss Integration Coefficients 


02715 24594 
06225 35239 
.09515 85117 
.12462 89713 
-14959 59888 
16915 65194 
18260 34150 
.18945 06105 
.18945 06105 
18260 34150 
16915 65194 
.14959 59888 
12462 89713 
09515 85117 
.06225 35239 
.02715 24594 


[3] for the interval (—1,1)). After using them the result must be multiplied by 
half the length of the interval in each case. 

Method of Checking. The tables were checked as follows: 

Table 1 was checked by using it to calculate 


f xJo(x)dx 
ri-1 


over each half-cycle and comparing the result with the known value of the integral 
[xJi (x) Ti, = xJi(xs) — xi-1Ji(xi-1), 


the values of J:(x,) being obtained from [5 ]. In each case agreement was obtained 
to within the accuracy warranted by the number of significant figures in the table. 
Table 2 was checked by using it to calculate 


_ Ji(x)dx = [—Jo(x) Bi. = Jo(yi-1) — Jol) 


over each half-cycle, where the y; are the zeros of J;(x). The results were calcu- 
lated to ten places of decimals and in each case were accurate to within two units 
in the last place. 
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TABLE 1 
x Jo(x) x Jo(x) 
.00000 00000 5.52007 81103 
1 .01274 44512 .99995 93952 1 5.53668 49893 .00564 19942 
. 2 .06664 37005 .99888 99624 2 5.60691 93443 .02928 21569 
3 .16156 67593 .99348 46849 3 5.73061 04883 .06978 62151 
. 4 .29410 48650 .97849 22054 4 5.90331 65740 -12301 12685 
5 45947 04869 .94791 40344 5 6.11879 91260 .18206 24818 
6 .65168 75524 .89661 10444 6 6.36927 09191 .23760 14372 
7 .86380 90708 .82197 96824 7 6.64567 95559 .27944 55482 
8 1.08816 85230 .72517 29892 8 6.93803 50453 .29921 03673 
9 1.31665 70348 .61136 06080 9 7.23577 09779 .29294 41222 
10 1.54101 64869 48883 00670 10 7.52812 64673 .26243 80011 
11 1.75313 80053 .36721 18144 11 7.80453 51041 -21449 26872 
12 1.94535 50709 .25549 12057 12 8.05500 68972 -15853 58154 
13. 2.11072 06927 .16052 05519 13 8.27048 94492 -10381 91971 
14 2.24325 87984 .08645 02192 14 8.44319 55349 .05742 97048 
15 2.33818 18573 .03506 03860 15 8.56688 66789 .02366 25625 
16 2.39208 11065 .00663 36641 16 8.63712 10339 .00451 20991 
2.40482 55577 8.65372 79129 
1 2.42133 49398 —.00854 11306 1 8.67035 68206 —.00450 94160 
2 2.49115 69581 —.04397 83230 2 8.74068 43400 —.02345 71287 
3 2.61412 19276 —.10342 59686 3 8.86453 95579 —.05612 15388 
4 2.78581 40782 —.17919 07027 4 9.03747 47347 —.09944 20369 
y 5 3.00003 15602 —.26006 26555 § 9.25324 31200 —.14809 16553 
6 3.24903 28661 —.33251 74168 6 9.50404 71589 —.19457 96817 
7 3.52381 87438 —.38328 09269 7 9.78082 24462 —.23045 12070 
8 3.81445 78522 —.40269 93971 8 10.07356 57393 —.24845 57197 
9 4.11044 58158 —.38757 11113 9 10.37169 66127 —.24484 61857 
10 4.40108 49242 —.34203 66196 10 10.66443 99058 —.22066 46191 
11 4.67587 08019 —.27600 29383 11 10.94121 51931 —.18130 82267 
al 12 4.92487 21078 —.20186 70491 12 11.19201 92320 —.13461 51679 
13 5.13908 95898 —.13110 13628 13 11.40778 76173 —.08848 00142 
14 5.31078 17404 —.07207 25099 14 11.58072 27941 —.04908 31439 
sd 15 5.43374 67099 —.02957 09059 15 11.70457 80120 —.02026 32192 
" 16 5.50356 87282 —.00562 57234 16 11.77490 55314 —.00386 81029 
j 5.52007 81103 11.79153 44391 
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x 


11.79153 44391 
11.80817 17029 
11.87853 45618 
12.00245 20168 
12.17547 40932 
12.39135 09018 
12.64228 09694 
12.91919 53359 
13.21208 57321 
13.51036 64157 
13.80325 68119 
14.08017 11784 
14.33110 12460 
14.54697 80546 
14.72000 01310 
14.84391 75860 
14.91428 04449 
14.93091 77086 
14.94755 90159 
15.01793 89756 
15.14188 65473 
15.31495 06747 
15.53087 99495 
15.78187 10027 
16.05885 26700 
16.35181 42497 
16.65016 74269 
16.94312 90066 
17.22011 06739 
17.47110 17271 
17.68703 10019 
17.86009 51293 
17.98404 27010 
18.05442 26607 
18.07106 39679 


TABLE 1—Continued 
Jo(x) 


.00386 
.02012 
.04823 
.08568 
.12799 
16873 
-20055 
.21699 
.21458 
.19401 
15987 
-11900 
.07837 
04355 
.01799 
00343 


— .00343 
— .01789 
— .04294 
— .07640 
— .11433 
— .15103 
—.17990 
— .19507 
— .19331 
—.17514 
-14458 
—.10779 
— .07109 
— .03954 
— 01635 
— .00312 


45915 
43811 
70070 
63326 
26729 
86824 
59550 
70655 
33409 
42035 
41093 
00109 
94720 
03186 
93958 
81248 


51301 
91458 
96191 
67698 
61880 
88400 
39278 
68842 
88429 
01558 
58966 
31239 
30831 
30642 
51985 
53477 
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x 


18.07106 
18.08770 
18.15809 
18.28206 
18.45514 
18.67110 
18.92213 
19.19915 
19.49215 
19.79054 
20.08354 
20.36056 
20.61159 
20.82755 
21.00063 
21.12460 
21.19499 
21.21163 
21.22828 
21.29867 
21.42265 
21.59575 
21.81173 
22.06277 
22.33981 
22.63284 
22.93126 
23.22428 
23.50133 
23.75237 
23.96835 
24.14145 
24.26543 
24.33582 
24.35247 


39679 
75348 
70515 
14541 
90820 
76781 
28135 
20923 
34535 
71443 
85055 
77843 
29197 
15158 
91437 
35463 
30630 
66299 
15866 
69809 
17345 
38153 
04440 
65401 
89501 
47770 
33838 
92107 
16207 
77168 
43455 
64263 


11799 * 


65742 
15308 
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Jo(x) 


-00312 
.01627 
.03908 
.06960 
.10428 
13794 
.16454 
17868 
17733 
-16087 
13297 
.09924 
.06551 
.03646 
.01509 
.00288 


— .00288 
— .01503 
— .03611 
— .06434 
— .09648 
—.12775 
— .15254 
— .16583 
— .16474 
— .14961 
— .12378 
— .09245 
— .06107 
— .03401 
— .01408 
— .00269 


28948 
88282 
92889 
69216 
43725 
52734 
42527 
44763 
10142 
71698 
88950 
99767 
96366 
96948 
18172 
47612 


27393 
12702 
16663 
88934 
81035 
47593 
57750 
05446 
72215 
08067 
03876 
99613 
91462 
63997 
19757 
23112 
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x 


24.35247 
24.36911 
24.43951 
24.56349 
24.73660 
24.95259 
25.20365 
25.48071 
25.77375 
26.07219 
26.36523 
26.64229 
26.89335 
27.10934 
27.28245 
27.40643 
27.47683 
27.49347 
27.51012 
27.58052 
27.70451 
27.87763 
28.09362 
28.34469 
28.62176 
28.91481 
29.21326 
29.50631 
29.78338 
30.03445 
30.25045 
30.42357 
30.54755 
30.61795 
30.63460 


15308 
74027 
66677 
82381 
98370 
83413 
82413 
58846 
78237 
28391 
47782 
24215 
23215 
08258 
24247 
39951 
32601 
91320 
56384 
75870 
38837 
20816 
88193 
82897 
64945 
96043 
59961 
91059 
73107 
67811 
35188 
17167 
80134 
99620 
64684 


TABLE 1—Continued 


Jo(x) 


00269 06195 
.01403 25020 
.03372 47328 
.06012 63461 
.09021 33604 
-11953 24224 
-14283 84354 
-15540 20499 
15450 96015 
.14042 06959 
-11625 87219 
.08689 57273 
05743 36567 
.03199 94405 
.01325 09005 
00253 38391 


00253 24019 
01320 95154 
.03175 59526 
05663 87732 
-08502 21095 
11271 68644 
13477 53040 
14672 13774 
-14596 96665 
.13273 89002 
-10995 96619 
.08222 81150 
.05437 12708 
03030 32039 
.01255 14239 
.00240 04013 
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x 


30.63460 
30.65125 
30.72165 
30.84564 
31.01876 
31.23477 
31.48584 
31.76292 
32.05598 
32.35444 
32.64750 
32.92457 
33.17565 
33.39165 
33.56477 


64684 
34327 
73180 
70254 
99856 
26651 
90422 
48688 
60401 
06419 
18132 
76398 
40169 
66964 
96566 


33.68876 93640 


33.75917 
33.77582 
33.79246 
33.86287 
33.98686 
34.15999 
34.37599 
34.62708 
34.90416 
35.19722 
35.49568 
35.78875 
36.06583 
36.31691 
36.53292 
36.70605 
36.83004 
36.90045 
36.91709 


32493 
02136 
75193 
28479 
50970 
16062 
87138 
02379 
17446 
89237 
96437 
68228 
83295 
98536 
69612 
34704 
57195 
10481 
83537 


Jo(x) 


00239 91633 
01251 61523 
03009 59796 
05369 50891 
08063 46686 
10694 79499 
12793 92673 
-13934 93292 
13870 48969 
12619 35228 
-10458 43845 
.07823 96934 
05175 15312 
02885 08462 
01195 21353 
.00228 60352 


.00228 49556 
.01192 16140 
.02867 
05116 72201 
.07686 29839 
10198 26451 
12204 77212 
13298 71140 
.13242 66679 
12052 95168 
— .09992 72078 
— .07478 03172 
— .04947 71799 
— .02758 90427 
—.01143 12041 
— .00218 65932 


17157 
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COND UF WD 


x 


36.91709 83537 
36.93374 59205 
37.00415 23534 
37.12814 65473 
37.30127 57720 
37.51728 62677 
37.76837 17300 
38.04545 75826 
38.33852 93584 
38.63699 47599 
38.93006 65358 
39.20715 23884 
39.45823 78507 
39.67424 83464 
39.84737 75711 
39.97137 17650 
40.04177 81979 
40.05842 57646 
40.07507 35355 
40.14548 08322 
40.26947 65473 
40.44260 78959 
40.65862 10416 
40.90970 95842 
41.18679 88362 
41.47987 42074 
41.77834 32704 
42.07141 86416 
42.34850 78936 
42.59959 64362 
42.81560 95819 
42.98874 09305 
43.11273 66456 
43.18314 39423 
43.19979 17132 


TABLE 1—Continued 


Jo(x) 


00218 56417 
01140 44585 
.02743 21950 
04896 56861 
.07357 53320 
.09765 02107 
.11690 15029 
12742 34635 
.12693 02235 
11556 53707 
.09584 13189 
.07174 25318 


.04747 84574 
.02647 94775 


.01097 29247 


-00209 90896 


.00209 82432 
.01094 92382 


.02634 06534 


04702 58304 
.07067 63342 
.09382 67282 
11235 56730 
12250 42137 
12206 57451 
11116 79231 
.09221 87488 
.06904 71411 


— .04570 38541 


.02549 38175 


—.01056 56711 


.00202 13118 


CONAN WD 


x 


43.19979 
43.21643 
43.28684 
43.41084 
43.58397 
43.79999 
44.05108 
44.32817 
44.62125 
44.91972 
45.21280 
45.48989 
45.74098 
45.95700 
46.13013 
46.25413 
46.32454 
46.34118 
46.35783 
46.42824 
46.55224 
46.72537 
46.94139 
47.19248 
47.46958 
47.76266 
48.06113 
48.35421 
48.63131 
48.88240 
49.09842 
49.27155 
49.39555 
49.46596 


17132 
96470 
76319 
45592 

76004 
28578 
38550 
58158 
40521 

60329 
42692 

62300 
72272 

24846 

55258 

24531 

04380 
83717 
64373 

49795 
28884 
73001 

42674 
72523 
14065 
19629 
63063 
68627 
10169 
40018 
09691 
53808 
32897 
18319 


49.48260 98974 


Jo(x) 


00202 
.01054 
.02536 
04529 
-06809 
.09041 
.10830 
-11811 
11772 
10723 
.08897 
.06663 
.04411 
.02461 
-01020 
.00195 


— .00195 
— .01018 
— .02449 
— .04375 
— .06577 
— .08735 
— .10465 
— .11416 
— .11380 
— .10369 
.08605 
— .06445 
— .04267 
— .02381 
— .00987 
— .00188 


05537 
45035 
98152 
96158 
50035 
97790 
19105 
38815 
07481 
69760 
80881 
43209 
43990 
06003 
06291 
15847 


08992 
15645 
89606 
05072 
72703 
88503 
74117 
40209 
89156 
54948 
66385 
79246 
99914 
32307 
09780 
86072 





on Ao FS WD 





x 

49.48260 

37 1 49.49925 
35 2 49.56966 
52 3 49.69366 
58 4 49.86680 
35 5 50.08281 
90 6 50.33391 
0S 7 50.61101 
15 8 50.90409 
81 9 51.20256 
60 10 51.49565 
81 11 51.77274 
09 12 52.02384 
90 13 52.23986 
03 14 52.41299 
91 15 52.53699 
A7 16 52.60740 
52.62405 

92 1 52.64070 
45 2 52.71110 
06 3 52.83510 
72 4 53.00824 
03 5 53.22426 
03 6 53.47536 
17 7 53.75245 
09 8 54.04554 
56 9 54.34402 
48 10 54.63710 
85 11 54.91420 
46 12 55.16529 
14 13 55.38131 
07 14 55.55445 
80 15 55.67845 
72 16 55.74886 
55.76551 








98974 
80712 
70711 
57858 
13228 
96942 
43111 
02663 
27276 
90110 
14723 
74275 
20444 
04158 
59528 
46675 
36674 
18411 
01048 
94851 
88697 
53420 
48803 
08536 
83057 
23502 
02460 
42905 
17426 
77159 
72542 
37265 
31111 
24914 
07550 


TABLE 1—Continued 
Jo(x) 


00188 
.00985 
.02371 
.04235 
06368 
08458 
10135 
11058 
-11026 
10048 
08340 
-06248 
04137 
.02308 
.00957 
.00183 


— .00183 
— .00955 
— .02299 
— .04107 
— .06177 
— .08206 
— .09835 
— .10732 
— .10702 
— .09755 
— .08098 
— .06067 
— .04018 
— 02242 
— .00929 
— .00177 


79842 
36907 
20292 
01575 
11286 
90824 
76521 
55944 
27587 
31520 
52240 
17090 
69630 
86479 
13465 
13567 


07874 
55769 
63548 
61891 
33832 
70662 
15046 
38205 
86465 
19386 
46590 
67598 
64126 
64133 
74390 
90152 


Conan kt OND 


x 


55.76551 
55.78215 
55.85256 
55.97656 
56.14970 
56.36572 
56.61682 
56.89392 
57.18700 
57.48548 
57.77857 
58.05567 
58.30676 
58.52278 
58.69592 
58.81992 
58.89033 
58.90698 
58.92363 
58.99404 
59.11804 
59.29118 
59.50720 
59.75830 
60.03539 
60.32848 
60.62696 
60.92005 
61.19715 
61.44825 
61.66427 
61.83741 
61.96141 
62.03182 
62.04846 


07550 
90943 
87941 
87414 
59995 
65182 
36311 
23408 
77154 
69658 
23404 
10501 
81630 
86817 
59398 
58871 
55869 
39261 
23295 
23003 
27250 
06495 
19998 
00793 
98556 
63585 
67579 
32608 
30371 
11166 
24669 
03914 
08161 
07869 
91902 


TABLES FOR EVALUATION OF INTEGRALS INVOLVING BESSEL FUNCTIONS 


Jo(x) 


.00177 
00928 
.02234 
03991 
.06002 
.07975 
09559 
-10433 
-10406 
09486 
.07876 
05901 
03909 
02181 


.00173 


— .00173 
— .00903 
— .02174 
— .03883 
— .05842 
— .07763 
— .09306 
— .10158 
— .10133 
— .09238 
— .07671 
— .05749 
— .03808 
— .02125 
— .00881 
— .00168 




















173 


84922 
29778 
17954 
06704 
74072 
79441 
78634 
45986 
33501 
31094 
32823 
96906 
30391 
80708 
57744 
09189 


04364 
24503 
01207 
90402 
15607 
33587 
32526 
19737 
15847 
50042 
51960 
13542 
43044 
66871 
34966 
65234 





174 TABLES 


COonaunr wd = 


TABLE 2 
x Ji (x) 


-00000 00000 
.02030 62503 .01015 26018 
.10618 61074 .05301 82575 
.25743 08620 12765 21143 
46860 91943 .22793 16731 
.73209 29378 .34206 47543 
1.03836 01743 45228 19201 
1.37634 19817 .53758 52259 
1.73382 29846 .57949 53889 
2.09788 29856 .56851 14490 
2.45536 39885 .50780 82279 10 
2.79334 57959 41190 88809 11 
3.09961 30324 .30107 17052 12 
3.36309 67759 19458 39510 13 
3.57427 51082 .10623 42640 14 
3.72551 98628 04329 54722 15 
3.81139 97199 .00819 97562 16 
3.83170 59702 
3.84857 90494 —.00678 05723 1 
3.91993 92273 —.03509 12559 2 
4.04561 30809 —.08323 73938 3 
4.22108 75619 —.14584 38197 4 
5 
6 
7 


Conant WH 


4.44002 41954 —.21441 93722 
4.69451 09307 —.27792 36119 
4.97535 02648 —.32471 41295 
5.27239 20719 —.34554 21188 8 
5.57490 05681 —.33641 08609 9 
5.87194 23752 —.29986 78712 10 
6.15278 17093 —.24400 22677 11 
6.40726 84446 —.17965 91930 12 
6.62620 50781 —.11727 67942 13 
6.80167 95591 —.06471 02564 14 
6.92735 34127 —.02661 42650 15 
6.99871 35906 —.00506 97404 16 
7.01558 66698 


x 


7.01558 
7.03232 
7.10309 
7.22774 
7.40178 
7.61893 
7.87134 
8.14989 
8.44450 
8.74454 
9.03916 
9.31770 
9.57011 
9.78726 
9.96130 
10.08595 
10.15673 
10.17346 
10.19016 
10.26076 
10.38511 
10.55873 
10.77535 
11.02715 
11.30502 
11.59892 
11.89823 
12.19213 
12.47000 
12.72180 
12.93842 
13.11204 
13.23639 
13.30699 
13.32369 


66698 
19654 
94235 
70376 
86125 
74325 
60560 
20857 
82807 
65243 
27193 
87490 
73725 
61925 
77674 
53815 
28396 
81351 
28485 
86758 
40273 
35623 
58100 
23585 
29316 
47008 
53706 
71398 
77129 
42614 
65091 
60441 
13956 
72229 
19363 


FOR EVALUATION OF INTEGRALS INVOLVING BESSEL FUNCTIONS 


Ji (x) 


00501 
.02606 
06226 
-11009 
-16354 
.21430 
-25311 
.27216 
.26752 
-24053 
19721 
.14615 
.09591 
.05314 
.02191 
.00418 


— .00416 
— .02167 
— .05192 
— .09212 
— .13743 
— .18092 
— .21471 
— .23197 
— .22906 
— .20682 
—.17021 
— .12656 
— .08328 
— .04624 
— .01910 
— .00364 


63154 
88611 
85746 
66693 
80168 
84782 
94883 
07618 
21198 
22625 
25700 
17084 
29504 
08898 
91616 
21325 


51430 
90724 
08365 
85148 
55137 
68681 
99134 
38712 
10599 
24276 
73009 
04959 
24257 
08298 
14741 
75651 








—=— i i ee 


ee ee a, le 
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TABLE 2—Continued 


3 Ji (x) x Ji(x) 
13.32369 19363 19.61585 85105 
154 1 13.34036 92371 .00363 92005 1 19.63252 14378 .00299 89710 
611 2 13.41090 14225 .01895. 70437 2 19.70299 28337 .01563 52051 
7146 3 13.53511 70816 .04546 56026 3 19.82710 14354 03755 37294 
693 4 13.70855 55312 .08082 85751 4 20.00039 04040 .06689 70210 
168 5 13.92495 18418 -12085 59397 5 20.21660 02095 -10026 96775 
1782 6 14.17648 57662 .15950 88072 6 20.46791 73449 -13270 31954 
1883 7 14.45406 65200 .18981 34949 7 20.74525 88608 -15837 97464 
1618 8 14.74766 17494 .20562 58827 8 21.03860 10500 -17208 93345 
1198 9 15.04666 02378 .20358 27840 9 21.33734 18412 .17088 33480 
1625 10 15.34025 54672 .18427 60601 10 21.63068 40304 15511 21291 
5700 11 15.61783 62210 .15200 47172 11 21.90802 55463 .12827 75447 
7084 12 15.86937 01454 .11324 27689 12 22.15934 26817 .09578 30326 
9504 13 16.08576 64560 .07464 19546 13. 22.37555 24872 .06325 40801 
3898 14 16.25920 49056 .04149 70116 14 22.54884 14558 .03521 86914 
1616 15 16.38342 05647 .01715 74777 15 22.67295 00575 .01457 70470 
1325 16 16.45395 27501 .00327 80206 16 22.74342 14534 .00278 66763 
16.47063 00509 22.76008 43806 
1430 1 16.48729 82913 —.00327 29251 1 22.77674 39258 —.00278 40720 
0724 2 16.55779 21583  —.01705 76459 2 22.84720 10186 —.01451 84126 
8365 3 16.68194 03340 —.04094 59048 3 22.97128 44308 —.03488 63684 
5148 4 16.85528 45588 —.07288 03022 4 23.14453 82281 —.06218 23121 
5137 5 17.07156 33065 —.10912 92847 5 23.36070 41509 —.09327 05856 
8681 6 17.32296 05786 —.14426 61835 6 23.61197 02781 —.12354 19341 
9134 7 17.60039 05295 —.17197 26270 7 23.88925 55038 —.14757 67318 
8712 8 17.89382 62558 —.18662 79783 8 24.18253 81552 —.16049 77977 
0599 9 18.19266 23056  —.18509 40367 9 24.48121 83130 —.15951 75470 
4276 10 18.48609 80319 —.16781 59539 10 24.77450 09644 —.14492 09337 
3009 11 18.76352 79828 —.13863 52301 11 25.05178 61901 —.11994 54048 
4959 12 19.01492 52549 —.10341 91552 12 25.30305 23173 —.08962 51276 
4257 13 19.23120 40026 —.06824 25400 13. 25.51921 82401 —.05922 29479 
8208 14 19.40454 82274 —.03797 23479 14 25.69247 20374 —.03298 98365 
ATA1 15 19.52869 64031 —.01570 98187 15 25.81655 54496 —.01365 90977 
5651 16 19.59919 02701 —.00300 24768 16 25.88701 25424 —.00261 16864 


19.61585 85105 25.90367 20876 











176 TABLES 


Con aAun Fr WH 


x 


25.90367 
25.92032 
25.99077 
26.11484 
26.28807 
26.50420 
26.75544 
27.03268 
27.32593 
27.62456 
27.91781 
28.19505 
28.44629 
28.66242 
28.83565 
28.95972 
29.03017 
29.04682 
29.06348 
29.13392 
29.25797 
29.43119 
29.64730 
29.89851 
30.17573 
30.46894 
30.76755 
31.06077 
31.33799 
31.58920 
31.80531 
31.97852 
32.10258 
32.17302 
32.18967 





FOR EVALUATION OF INTEGRALS INVOLVING BESSEL FUNCTIONS 


20876 
93474 
67744 
31639 
31930 
94607 
11174 
83032 
07199 
99027 
23194 
95052 
11619 
74296 
74587 
38482 
12752 
85349 
41780 
47675 
91152 
23308 
76205 
48930 
51694 
91244 
93216 
32766 
35530 
08255 
61152 
93308 
36785 
42680 
99110 


TABLE 2—Continued 
Ji (x) 


00260 
.01361 
.03271 
.05834 
.08755 
11605 
13872 
15097 
15015 
13650 
11305 
08452 
05587 
.03113 
01289 
.00246 


~- .00246 
— .01285 


— .03090 


— 05513 


— .08278 


— .10977 


— .13129 


— .14296 
— .14227 
— .12941 
— .10722 
— .08020 
— .05304 
— .02956 
— .01224 
— .00234 


96493 
13779 
75400 
28100 
98843 
08696 
26066 
45446 
74510 
89575 
33742 
19696 
69690 
74830 
55290 
60528 


43999 
56731 
89818 
72476 
50383 
68854 
35471 
90393 
43786 
16541 
87815 
27911 
14604 
61656 
73498 
23778 


Conan r WN 


x 


32.18967 99110 
32.20633 43684 
32.27676 99432 
32.40081 54597 
32.57401 63444 
32.79011 62490 
33.04130 56382 
33.31850 61796 
33.61169 92609 
33.91028 82003 
34.20348 12816 
34.48068 18230 
34.73187 12122 
34.94797 11168 
35.12117 20015 
35.24521 75180 
35.31565 30928 
35.33230 75501 
35.34896 11121 
35.41939 29002 
35.54343 17477 
35.71662 33207 
35.93271 16074 
36.18388 74921 
36.46107 31306 
36.75425 04492 
37.05282 33358 
37.34600 ~06544 
37.62318 62929 
37.87436 21776 
38.09045 04643 
38.26364 20373 
38.38768 08848‘ 
38.45811 26729 
38.47476 62348 


Ji(x) 


00234 
.01221 
.02937 
.05240 
.07871 
.10442 
-12494 
13611 
13551 
12331 
.10222 
.07648 
.05059 
.02821 
.01168 
.00223 


— .00223 
— .01165 
— .02804 
— .05004 
— .07519 
— .09978 
— .11943 
— .13016 
— .12963 
—.11801 
.09785 
— .07324 
— .04846 
— .02702 
— .01119 
— .00214 


09996 
33904 
08010 
83042 
53604 
22839 
39101 
63937 
63765 
84545 
15009 
52251 
84188 
12428 
81140 
56390 


44653 
87150 
16792 
84236 
26901 
23087 
51446 
38763 
88064 
27490 
65260 
13135 
48482 
71631 
91603 
22896 
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TABLE 2—Continued 


x Ji(x) x Ji(x) 
38.47476 62348 44.75931 89977 
996 1 38.49141 91040 .00214 12735 1 44.77597 08809 .00198 52730 
904 2 38.56184 79626 .01117 33773 2 44.84639 55686 .01036 06792 
010 3 38.68588 16508 .02687 81611 3 44.97042 19116 .02492 87731 
042 4 38.85906 60202 .04798 12388 4 45.14359 60250 .04451 51717 
604 5 39.07514 53189 .07210 44757 5 45.35966 25275 .06692 15430 
839 6 39.32631 07563 .09571 08701 6 45.61081 30909 .08887 03803 
101 7 39.60348 48656 .11459 65682 7 45.88797 07859 -10645 75916 
937 8 39.89664 99899 .12493 01313 8 46.18111 85490 -11611 55971 
765 9 40.19521 04577 -12446 56086 9 4647966 13360 -11574 22770 
545 10 40.48837 55820 -11333 82116 10 46.77280 90991 -10544 63168 
009 11 40.76554 96913 .09400 72264 1i 47.04996 67941 .08750 13044 
251 12 41.01671 51287 .07037 82521 12 47.30111 73575 .06553 43925 
188 13 41.23279 44274 -04658 04345 13 47.51718 38600 .04338 95839 
428 14 41.40597 87968 .02598 07709 14 47.69035 79734 .02420 77320 
140 15 41.53001 24850 .01076 68903 15 47.81438 43164 .01003 40865 
390 16 41.60044 13436 .00205 97434 16 47.88480 90041 .00191 97690 
41.61709 42128 47.90146 08872 
653 1 41.63374 65353 —.00205 88518 1 47.91811 24121 —.00191 90607 
150 2 41.70417 30810 —.01074 40331 2 47.98853 55848 —.01001 56479 
792 3 41.82820 26960 —.02584 84372 3 48.11255 92598  —.02410 07718 
236 4 42.00138 13780 —.04615 08277 4 48.28572 96479 —.04304 20054 
901 5 42.21745 35808 —.06936 81412 5 48.50179 15025  —.06471 68956 
087 6 42.46861 07699 —.09210 05548 6 48.75293 66631 —.08595 79773 
446 7 42.74577 57769 —.11030 23920 7 49.03008 83960 —.10298 88602 
763 8 43.03893 12737 —.12028 11784 8 49.32322 98530 —.11235 50074 
064 9 43.33748 19369 —.11986 63996 9 49.62176 62178 —.11201 66719 
490 10 43.63063 74337 —.10917 88005 10 49.91490 76748 —.10207 23609 
260 11 43.90780 24407 —.09057 93974 11 50.19205 94077 —.08471 72272 
135 12 44.15895 96298 —.06782 68662 12 50.44320 45683 —.06345 97884 
482 13 44.37503 18326 —.04490 01454 13 50.65926 64229 —.04202 19605 
631 14 44.54821 05146 —.02504 72789 14 50.83243 68110 —.02344 73534 
603 15 44.67224 01296 —.01038 11283 15 50.95646 04860 —.00971 96893 
896 16 44.74266 66753 —.00198 60641 16 51.02688 36587 —.00185 97015 


44.75931 89977 51.04353 51836 











178 TABLES 


Conant WH 


x 


51.04353 
51.06018 
51.13060 
51.25462 
51.42779 
51.64385 
51.89499 
52.17214 
52.46527 
52.76380 
53.05694 
53.33409 
53.58523 
53.80129 
53.97445 
54.09848 
54.16890 
54.18555 
54.20220 
54.27262 
54.39664 
54.56980 
54.78586 
55.03700 
55.31414 
55.60727 
55.90580 
56.19893 
56.47607 
56.72721 
56.94326 
57.11643 
57.24045 
57.31087 
57.32752 





FOR EVALUATION OF INTEGRALS INVOLVING BESSEL FUNCTIONS 


51836 
64126 
83339 
98048 
71154 
51303 
58277 
26353 
88827 
99421 
61895 
29971 
36945 
17094 
90200 
04909 
24122 
36411 
46228 
54982 
51273 
98663 
46726 
16405 
43323 
62264 
28526 
47467 
74385 
44064 
92127 
39517 
35808 
44562 
54379 


TABLE 2—Continued 


Ji (x) 


-00185 90619 
00970 29509 
.02335 01746 
.04170 60842 
.06271 67691 
08331 44272 
.09983 86147 
10893 77477 
.10862 92546 
-09900 29593 
.08218 31360 
06157 05817 
.04077 60577 
.02275 44296 
.00943 31183 
.00180 49432 


— .00180 43618 
— .00941 78326 
— .02266 56242 
— .04048 73357 
— .06089 13712 
— .08090 07468 
— .09696 09091 
—.10581 45509 
—.10553 17575 
— .09619 48735 
— .07986 37153 
— .05984 06981 
— .03963 48225 
— .02211 95390 
—.00917 04951 
—.00175 47553 


Conant ON 


x 


57.32752 
57.34417 
57.41459 
57.53861 
57.71177 
57.92782 
58.17896 
58.45610 
58.74923 
59.04775 
59.34088 
59.61802 
59.86915 
60.08520 
60.25836 
60.38238 
60.45280 
60.46945 
60.48610 
60.55652 
60.68054 
60.85370 
61.06975 
61.32088 
61.59802 
61.89114 
62.18966 
62.48279. 
62.75992 
63.01105 
63.22710 
63.40027 
63.52428 
63.59470 
63.61135 


54379 
62108 
62035 
42778 
68459 
89436 
27631 
19805 
01998 
30834 
13027 
05201 
43396 
64373 
90054 
70797 
70724 
78453 
84404 
76810 
44308 
51497 
49401 
60776 
23352 
74241 
71197 
22086 
84662 
96037 
93941 
01130 
68628 
61034 
66985 


Ji (x) 


00175 42236 
00915 64627 
.02203 79680 
03936 95712 
.05921 66524 
.07868 54349 
.09431 85737 
-10294 55141 
-10268 50402 
-09361 30005 
.07773 02633 
.05824 89261 
03858 43703 
.02153 50060 
00892 86590 
.00170 85351 


—.00170 80464 
— .00891 57171 
— .02145 97360 
— .03833 95708 
— .05767 29573 
— .07664 27230 
— .09188 11714 
—.10029 79275 
—.10005 69840 
— .09122 85052 
—.07575 91800 
— .05677 78096 
— .03761 32667 
— .02099 45050 
— .00870 50043 
—.00166 57858 
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Example of the Use of the Tables. The example 


f Jo(x)dx = 1 
0 


given in [1] is repeated here partly since an error in Watson ([6], p. 752) is 


2236 revealed, and partly as an extension of the tabulation of 
1627 ce 
a f Jo(x)dx 
) mm 
9524 F r 
1349 to a larger value of m than that obtainable from Watson ([6], p. 752). 
5737 Applying the integration coefficients directly to the values of Jo(x) in Table 1 
5141 or 
)402 Zn 
005 nf” Jo(x)dx a? at at as 
2633 1 1.47030 0043 
9261 
3703 2 .80145 4213 
060 3 .59932 2516 
on 4 49904 9621 
5 .43653 5113 
~ 6 .39282 2560 
i360 7 .36005 6836 
en 8 .33432 0981 
7230 9 .31341 5072 
[714 10 .29599 6007 
9275 
11 .28119 1319 
ond — 12784 903 
5052 12 .26840 6416 15987 49 
1800 —11186 154 — 30760 6 
3096 13 .25722 0262 12911 43 76952 
— 9895 011 — 23065 4 —23112 
2667 14 24732 5251 10604 89 53840 7965 
050 — 8834 522 —17681 4 —15147 — 3058 
043 15 .23849 0729 8836 75 38693 4907 1263 
— 7950 847 — 13812 1 — 10240 —1795 —489 
858 16 .23053 9882 7455 54 28453 3112 774 
— 7205 293 — 10966 8 — 7128 —1021 
17 .22333 4589 6358 86 21325 2091 
— 6569 407 — 8834 3 — 5037 
18 .21676 5182 5475 43 16288 
— 6021 864 — 72055 
19 .21074 3318 4754 88 
— 5546 376 


20 .20519 6942 


Here x» denotes 0. 
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Now, choosing the first advancing row of differences we have (by Euler's 
transformation) approximately 


f Jo(x)dx = 1.47030 0043 — 0.80145 4213 + .59932 2516 — .49904 9621 
0 


+ 43653 5113 — .39282 2560 + .36005 6836 — .33432 0981 
+ .31341 5072 — .29599 6007 + (1/2)(0.28119 1319) 

+ (1/4) (0.01278 4903) + (1/8) (0.00159 8749) 

+ (1/16) (0.00030 7606) + (1/32) (0.00007 6952) 

+ (1/64) (0.00002 3112) + (1/128) (0.00000 7965) 

+ (1/256) (0.00000 3058) + (1/512)(0.00000 1263) 

+ (1/1024) (0.00000 0489) = 0.99999 9992. 


It should be noted that the first figure in the table 
(4) 1.47030 00434 = f * Jo(x)dx 
Jo 


disagrees with the figure for this integral quoted in [1], where it is taken from 
Watson ([6], page 752). This integral was therefore calculated independently 
by interpolation, 


f  Jolt)dt = { To(t)dt + hIo(x) -— = a Jie) += = ( 2 sl u(x) ), 


(Taylor’s theorem) using tables of Jo(x), Ji(x) [4] and a table [3] of 
So? Jo(t)dt. For this x was taken to be 2.4 and h to be 0.00482 55577, so that 
x + h = 2.40482 55577 is the first zero of Jo(x). The result obtained agreed with 
(4) to ten decimal places, and this confirmation reveals an error in Watson [6] 
where the value 


1 ( 
; f Jo(x)dx = 0.73522 08 


is quoted. This value should be corrected to 0.72515 00. 


I. M. LONGMAN 
Department of Applied Mathematics, 


The Weizmann Institute of Science, 
Rehovoth, Israel 


R JL. M. Loneman, ‘Note on a method for ome infinite integral§ of oscillatory func- 
tions,’ the Phil. Soc., Proc., v. 52, 1956, 64. 

2. T. J. I'A. BRomwicu, ‘An Introduction to ‘the Theory of Infinite Series, The Macmillan Co., 
New York, 1942. 

3. NBS Applied Mathematics Series, No. 37, Tables of Functions and of Zeros of Functions, 
U. S. Gov. Printing Office, Washington, D. es 1954. 

4. HARVARD UNIVERSITY, Computation Laboratory, Annals, v. 3: Tables of the Bessel Functions 
of the First Kind of Orders Zero and One, Harvard University Press, Cambridge, Mass., 1947. 

5. BRITISH ASSOCIATION For THE ADVANCEMENT OF SCIENCE, ‘Committee on Mathematical 
Tables, Bessel Functions, v. 6, ‘‘Zeros of Jo(x) and Ji(x),” Cambridge: at the University Press, 
1950. 

6. G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge: at the University 
Press, ee 

7. Davis & P. RaBinowitz, ‘“Abscissas and weights for Gaussian quadratures of high 
in NBS Jn. of Research, v. 56, Research Paper 2645, 1956, p. 35-37. 








’ 


rs 


om 
tly 


‘of 
hat 
ith 
6] 


inc- 


ess, 
sity 


igh 





TIME-STEP FOR COMPUTATION OF ORBITS BY NUMERICAL INTEGRATION 181 


On the Time-Step to be Used for the Computation 
of Orbits by Numerical Integration 


W. J. Eckert [1] has adapted Cowell’s method of numerical integration to the 
determination of orbits on punched card machines. Eckert, Brouwer, and 
Clemence [2] have used Cowell’s method on a large-scale computer. We shall 
assume that Cowell’s method as modified by Eckert is the method of numerical 
integration best suited for the determination of orbits on large-scale computers. 
What we are concerned with in this paper is how to make use of this method in 
such a way that we do the least work (have the fewest arithmetic operations) per 
unit of advance in the time. Our results apply also to the computation of orbits 
on ordinary desk computers, but additional factors, such as the work involved in 
transcribing, must be considered. 

In Cowell’s method as modified by Eckert the attractions and the central 
differences of the attractions which are required at the “‘new time-step” are 
obtained from the values at previous time-steps by an extrapolation process. 
This extrapolation process may be shown to be equivalent to extrapolating the 
attractions ahead in the time using Newton’s backward difference interpolation 
formula, so that Cowell’s integration formula may be written in terms of backward 
differences instead of central differences. Let x,', 1 = 1, 2, 3, be rectangular co- 
ordinates and let X,‘ be the attractions at time t = mAt. Then we obtain for x‘n4: 
in terms of backward differences, 


R? 


(1) wings = "X44, + AP ( = 


1 19 3 
X,§ + — VX, + — X,! + — WX, 
+n + VX. + Vx 


863 275 33,953 
4 Ay v5 nt : 
12,096 ve * 4,032 waa 518,400 


8,183 3,250,433 4,671 
; v7 A Bs, aceite a a V8X,,' : 
129,600 Xa" + 53,222,400 + 78,848 











os v8X, 


am 








v°x,' 





301,307 ,139,941 
vio 3° --+ } 
5,230,697 472,000 alltel ) 


(2) “Xi. = 2"X,§ — "X41 + APX,'. 
Here V is a backward difference operator, i.e., 
VX,' = X,' — X41, VX! = VX,' — VX'4-1, ete. 


The ”’X,,* defined by equation (1) are called ‘‘second-summations” [2]. In order 
to start the integration, "Xo and "X‘*_, and all the backward differences for 
n = 0 must be obtained in such a way that the initial conditions are satisfied. 
One way of doing this is to start the integration using a less accurate but “‘easy- 
starting” method with time-step a product of an appropriate negative power of 
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two times At. Then from this solution backward differences can be evaluated 
and then the second summations for starting Cowell’s method can be obtained 
from (1). 

Let m be the order of the last backward difference retained in (1). We shall 
call m — 2 the order of the integration formula (1) and we label the coefficient 
of the last backward difference retained, c,,. For a computer with large memory 
we may construct a sequence of instructions for integrating with (1) and (2) in 
which m and At are parameters, i.e., we may compute using any value of m 
between, say, one and one hundred and any value of At. We would like to choose 
the combination of m and At which minimizes the work performed by the com- 
puter in obtaining the orbit to a prescribed degree of accuracy. Assuming that 
most of the work is involved in computing the attractions and not in evaluating 
the backward difference series we would at first expect to choose At close to the 
limiting value for which the backward difference series in (1) converges. For a 
circular orbit this limiting value of At may be shown to be the time required to 
traverse an arc of 30° so that we would choose approximately 12 steps per period. 
However, as will be shown, with 12 steps per period, the integration would become 
unstable for m greater than 6, so that the radius of convergence of the backward 
difference series in (1) is not what reaily determines how large a At we should use. 
It turns out that if m is large enough or if At is large enough, spurious solutions 
of (1) and (2) initiated by rounding errors grow exponentially. In this case we 
say that the integration procedure is unstable. The spurious solutions are solutions 
of (1) and (2) which have no analogue in the differential equations of motion, i.e., 
which vanish as At — 0. Brouwer [3] has discussed the accumulation of rounding 
error when we integrate using (1) and (2) neglecting the backward difference 
series, i.e., when the spurious solutions damp out. Here we wish to estimate the 
value of At at which the integration first becomes uhstable as a function of m. 

Applying V? to equation (2) and substituting AfX,* for V?’’X+,,, we obtain 


(3) Xana = Qxtai — 24 + AP(XaS + VX, + --- + 6,.0"X,'). 


The integration formula (3) is derived by Collatz [4] with coefficients through 
the fifth order and is attributed by Collatz to Stérmer [5]. The solution of (1) 
and (2) depends on the starting values of ’’X,,‘ and 'X,,‘ as well as m — 1 values 
of the coordinates. The solution of (3) depends on m + 1 values of the coordinates. 
Hence the solution of equations (1) and (2) depends on the same number of 
arbitrary constants as the solution of (3). By choosing the starting conditions 
properly any solution of (1) and (2) is a solution of (3) and vice versa. 

Let us assume that we have a solution to (3) and that then we make a small 
perturbation in this solution. Let e,‘ be the difference between the coordinates for 
the new perturbed orbit and the old orbit. Let us suppose that the perturbation 
is so small that terms in (e,*)? are negligible. Let us assume further that we have 
only one body in the field of a fixed mass M. Then X,' = — yMx,‘/r,*. Let 
us use the notation X,,*? = 0X,‘/dx,’. We obtain equations for e,* by subtracting 
the equations for the perturbed orbit from the equations for the unperturbed 
orbit and then expanding the attractions for the perturbed orbit about the attrac- 
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tions for the unperturbed orbit keeping first order terms in «,‘ only. We obtain 


3 
(4) €'n+1 — ye + Cant = Af 2 LX nt tenf + TeV? (X,*4€,/) + --- 
7=1 
+ CmV™(X n**4€n/) }. 


The X,‘ are understood to be given functions of the coordinates x,‘ of the 
unperturbed orbit. Let the components of the orthonormal characteristic vectors 
of the 3 X 3 matrix X,‘/ belonging to characteristic values U,, V,, W, be u,‘,v,', 
w,' so that >>; X,*4u,i = U,u,', 5; Xn, = Vion‘, and >; X,*4w,i = W,w,'. 
Written out explicitly, the matrix X,‘/ appears as follows: 


( 1\2 1 2 1 om 7 
“iste ) CAR) eRe 
Tr Tr Tr Tn Tn 
M m a 2 2 at ne 
o weattla(e\(2) —r4a(2) 92) 
fs Tr Tn Tn Tn Tn 
3 1 3 2 3\2 
ead. eh eed 
{ Sa So. Tn Tn Tn 
By solving the determinantal equation belonging to this matrix we find that 


X,* has a single root 2y_M/r,? and a double root —yM/r,'. Label the single root 
U, and its corresponding characteristic vector u,*. We obtain for u,', 











(6) u,* - Sa*/ Par 


i.€., 4,‘ is a unit vector in the direction of the radius vector. The remaining charac- 
teristic vectors of X,,‘:’, v,* and w,‘, may be chosen arbitrarily in the plane per- 
pendicular to the radius vector. We expand e,‘ in u,', v,‘, w,*. Let the expansion 
coefficients be an, Ba, Yn SO that 


(7) en! > GnUn* + Badn* + YnWn*. 


In order to carry our analysis further we have found it necessary to make 
the assumption that the orbit is circular. Our results will hold approximately for 
non-circular orbits in terms of the local radius of curvature provided that the 
radius of curvature of the orbit does not change too rapidly. Let r, @ be polar 
coordinates in the plane of motion. Let @, be the value of @ corresponding to x,'. 


M ‘ 
It is easily shown that A@ = 7 Af, where A@ is the angle between two suc- 


cessive points on the orbit separated in time by At. We assume that the unper- 
turbed orbit lies in the x, y plane and that y = 0 and x = r when ¢ = 0. Then 


(8) x,) =x, =rcosndé, x,” = y, = rsin ndé, 
cos nAé@ 
(9) u, =| sin nA@ 


0 
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We choose for »,‘, w,', 


— sin nAé 0 
(10) v, = cosnA6 |, w,={| 0 
0 1 


Let S be the rotation matrix for rotating a vector through angle Aé@ in the 
(r, 0) plane. Then 


cos A@ —sinAéd 0 
(11) S={[sinA@ cosA@ 0 
0 0 1 


We have the following relations involving u,, v,, w,, and S: 


Un = S"to, Ua = SU, Wa = SW, 
UioS"Uy = ToS"Vo = cos nAO, 
(12) DoS"Uo = — MoS" = sin NAO, 
WoS*Uy = WoS"Vo = tioS"Wo = ToS" = O, 
WoS"Wo = 1. 


We obtain equations for the expansion coefficients an, Ba, yn by first substitut- 
ing (7)-into (4), then substituting for u,, S"uo; for v,, S"vo; and for wy, S*wo, 
and finally multiplying the equations obtained from the preceding steps suc- 
cessively by tipS—", toS-", WoS-". We then obtain 


(13) a@n41 cos Ad — 2a, + an_1 cos A@ — Bayi sin AO + B,_1 sin Ad 
- = 2A@A (an) + APB(B,), 


(14) anyi sin A? — ap_: sin A@ + 8,41 cos A@ — 28, + Bn_1 cos AO 
= 2APB(a,) — APA (B,), 


(15) Yn+1 — 2¥n + Yn—-1 = — A®Pc(yn), 


where A and B are the real and imaginary parts of J and 
(16) T (cn) = e7™L,eind 4+ DV? (wneind?) +... + 6,0" (wneimd?) T, 


(17) C(¥n) = Yn + PxV"¥n + + +> + CmV™Yn. 3 
In (13) and (14) the argument of A and B indicates that a, or 8, is to be sub- 
stituted for w, in [(w,). We consider first equation (15). Equation (15) has m 
solutions of the form y, = A", \ being obtained by substituting y, = A” in (15), 
multiplying by \"-*—', and solving the resulting polynomial of degree m in X. 
It may be shown that the two roots \ which approximate the solution of the 
variational differential equation obtained from (15) when A@— 0 lie almost on 
the unit circle close to +1, whereas if A@ is large enough some roots have negative 
real parts. As A@é is increased, one of these roots with negative real part emerges 
from the unit circle at \ = — 1 on the real axis. The value of A@ for which A = — 1 
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is the largest value of A@ for which the modulus of all roots is less than or equal to 
one, and is the largest value of A@ for which the integration is stable. Let A@,” 
be this largest value of A@ for which the integration formula of order m is stable. 
Substituting A = — 1 in (15) we obtain for Aé,”, 


2 3 


2 =~? 
(18) A6, w2[1t 7474... tearm] : 


Values of 46,” are tabulated as a function of m in Table 1 at the end of this paper- 

Next consider equations (13) and (14). These equations have solutions of the 
form a, = ad", B, = bd", X being obtained by substituting aX", bd" for a,, 8, in 
(13) and (14); multiplying the resulting equations by \”~*~", setting the deter- 
minant of coefficients of a, 6 to zero; and sclving the resulting determinantal 
equation of order 2m. We assume as before that \ = — 1 gives the largest value 
of Aé for which all roots have modulus less than one. A proof of this statement for 
the polynomial equation obtained from (13) and (14) is not easy, and, indeed, the 
statement may not even be true. However, it is certainly true that if \ < — 1, 
the integration is unstable, so that the assumption A = — 1 gives us a bound on 
Aé@ which may not be exceeded. We have had an opportunity to test the stability 
of various order time-step combinations on a computer and have found that the 
equations of motion do become unstable at values of Aé@ close to those predicted 
here, for the order time-step combinations tested. This gives us confidence that 
as A@ is increased, the first “spurious” root to leave the unit circle must emerge 


at a point close to minus one. Substituting w, = #(—1)" and making use of the 
relation 


46 \? 
V?(—1)eind? — (—1)neiln—vi2)a0 (2 cos 7 ) » 


we obtain 


(19) I(w) = (—1)* [1 + pre7 4? (2 cos + pee ii0e (2 cos “* wee 


2 
+ cae (2 cos Mf) |e. 


Then substituting a, = a(—1)", 8, = 6(—1)" in (13), (14), and multiplying 
through by (—1)* we obtain 


2 
(20) — 2(1 + cosAé)a = 2Af [: + 7s cos Ad (2 cos “*) 


; m 
+ dy cos 440 (2 cos “#) ee + ¢q.cos 40(2 cos“) Jo 


2 3 
— a6 | zr sin 20 (2 co “*) + drsin $40 (2.cos 2) a 


+ Cm» sin - Aé (2 cos a b, 
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Ae \? 
(21) — 2(1+ cos Aé)b = — AP] 1 + 4 cos AO 2 cos > 


r m 
+ #s cos 4.40 (2 cos “? ) ease + &m cos a9 (2 cos ) |. 


, Ae \? F Ae \3 
— 208 | asin 10 (2 cos“) + # sing (2 cos“) er 
Aé \™ 
+ en sin ™ 40 (2.008 2) Je. 


We solved (20) and (21) for Aé@ for several values of m. The results obtained 
are tabulated under the heading A@™,,, in Table 1. These solutions were obtained 
by successive approximation starting with A@,” and then substituting for the 
trigonometric functions in (20), (21); then solving the determinantal equation 
for Aé@, then resubstituting for the trigonometric functions, etc. Two or three 
iterations were sufficient to obtain A@”.,g to the number of places given in Table 1. 

We note that Aé™. is close to A@,”. We use Aé,” as the limiting value of Aé@ 
for stable integration in subsequent discussion. This approximation is justified 
since in practice we would not want to choose Aé close to its limiting value because 
this would make the rounding error large, and it also is justified because orbits 
obtained by numerical integration are rarely really circular so that we can only 

2a 
Ae,” 
N,, is the smallest number of steps per period giving stable integration of order m. 

In order to make use of our results we require an estimate of the truncation 
and rounding error. We assume that At has been chosen small enough so that the 
backward difference series in (3) converges. Let y,‘, Y,‘ be coordinates and attrac- 
tions corresponding to the exact solution of the differential system. Then y,', 
Y,‘ satisfy 


know the limiting value of A@ approximately. In Table 1 we tabulate N,, = 





inti — 2Yn' + Yat 





= ia crV* Y,,' => y a c.V* Y,,*. 


2 
At k=0 k=m+1 


(22) 


Our numerical solution x,‘ satisfies 





Saas ae 2xn' + Xing , 
—_ z CxV*Xn* = pa 


23 
(23) Ar a 


Here p, is the rounding error incurred at step ”. In order to keep as much accuracy 
as possible without performing extended accuracy multiplication, the term in 
parenthesis in (3) should be evaluated with such scaling as to give all precision 
possible using “‘single-precision” arithmetic, then extra digits should be carried 
in the summation of the coordinates in (3). Now let «,' = y,' — x,*. Assuming 
that ¢«,‘ is small, we obtain from (22) and (23) (in manner analogous to the 
derivation of (4)) 





ts i Re.t 4 of, = s aT ge 
Cnt én’ + e'n+1 d xy > CEV* (Vani den!) = — pat CupiV"t'Y,!. 


(24 
(24) At ae 
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In (24) we retain only the first term on the right-hand side of (22). ¢n4:V"*' Y,' 
is called the local truncation error. For the circular orbit defined by (8) 





M Aé@ m+1 
vt Y,! = rE (2 sin #) cos (wt, _ Ym); 
95 M A@ m+1 
(25) very? = — 1A (2 sin ) sin (wl, — ¢m), 
v=tiy3 = 0, 
1 Aé m+1 
where w = 4 px and ¢» equals — Z (40 — x). We let Anyi = Cnt (2 - z ) 


As can be seen by examining equations (22) and (23), if the problem is scaled so 
that (Y,*)max = yM/r occupies the full number size of the computer, then the 
local truncation error is the same order as the rounding error when A,,;; is the 
same order of magnitude as p,. This assumes that extra accuracy is carried in the 
coordinates as mentioned previously. A,,;: is tabulated in Table 1. 

We may now apply our results and choose an appropriate time-step size and 
order, m, for various particular cases. For example, suppose we wish to integrate 
for orbits of the planets on a computer with a 14 decimal-digit number size and 
we wish to obtain the orbits with all the precision of which the machine is capable 
without performing extended accuracy operations other than extended accuracy 
summations. We see from Table 1 that we are going to need between 60 and 70 
steps with m = 12 to keep the truncation error small enough and that then we 
are too close to the unstable value of A@, so that we must use more steps, say 100 
steps per period. But then the truncation error is sufficiently small with m = 10, 
11, or 12. m = 11 and 100 steps per period would appear to be a good compromise. 
In the integrations of the orbits of the five outer planets, Eckert, Brouwer, and 
Clemence [2] used m = 11 and they had approximately 113 steps per period in 
Jupiter, the planet with the smallest period. Hence this theoretical analysis gives 
results in good agreement with what has been found appropriate in practice. 

Many computers have a 10 decimal-digit word size. On these m = 10 with 60 
steps per revolution would appear to be a good choice. 

We see that we would not expect to use very high orders (m > 14) under any 
practical circumstances. An extreme case would be one in which we wished to 
follow a very large number of revolutions, say 10" revolutions. We would then 
wish to perform extended accuracy computations on, say, a 14 digit machine 
thus obtaining effectively 28 decimal-digit number size. In this case, to keep the 
local truncation error below 10-** and to be assured of stability we would want to 
choose about 180 time-steps per period and m = 13. Assuming we could perform 
the computations for one time-step in a millisecond (this is probably somewhat 
faster than we could expect from any computer today) it would require 56 years 
to complete the computation. Hence even in an extreme case m is less than 14. 

Of course we do not always wish to integrate obtaining the maximum precision 
of which a particular computer with its given number size is capable. Instead, we 
often wish the error in the orbit to be less than a given bound over the total 
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interval of the integration. In order to be assured of this it is necessary to estimate 
the accumulated truncation error and the accumulated rounding error. From our 
prescribed bound in the maximum allowable accumulated truncation error we 
obtain the maximum allowable At for each order m. Then for least work by the 
computer we choose the (m, At) combination which gives the largest value of At 
for which the integration is stable according to Table 1. We obtain the minimum 
allowable number size from our estimate of the accumulated rounding error. Then 
we have obtained the maximum Aft and minimum number size which may be used 
to solve our problem using Cowell’s method. 

Let us estimate the accumulated truncation error. It is easily seen that equa- 
tion (24) approximates the differential system 


3 
(26) @ei/dt? — 5 Vii = — p + CuyrV™¥! 
j=1 
where e*, Y*4, p, and ¢n4i1V" Y‘ now stand for continuous quantities. For the circu- 
lar orbit we can obtain exactly a particular solution of equation (26) corresponding 
to the local truncation error. It is 


(27) __e* (truncation) = rAniil(} CoS gm + 20 sin gn)u' — 34 sin g,v* ]. 


Here ¢» is the phase factor defined following equation (25) and u‘ and v‘ are de- 
fined in equations (9) and (10) except that here ”Aé is replaced by a continuous 
variable 6 = wt. We see that after a large number, P, of periods, the last term 
in (27) dominates and the amplitude of the truncation error approaches 
|6x? sin gmP?Am4i7|. The way the rounding error accumulates can also be in- 
vestigated starting with (26). One is led to the same result obtained by Brouwer 
[3], namely that after steps the root-mean-square rounding error increases by 
a factor of n!. The method used here to estimate the accumulated truncation error 
is similar to that discussed by Rademacher [6 ]. 


TABLE 1. Values of A0,™, AO"a,8, Nm, Amsi1 for m = 6(1)14 
m Ad,” AO”... 8 Nan Am+1 M% Ae,” Ad”... 8 Nun An+1 


6 .6252 .629 10.05 10% 11 .1260 49.86 10-” 

7 4593 13.68 10 12 .0905 .0893 69.39 10-% 

8 .3346 .3211 18.78 10-* 13.0650 96.70 10-8 

9 .2424 25.92 10° 14 .0466 134.88 10-*! 
10 .1750 .1701 35.90 10-” 


Note: Values listed for m = 13, 14, and A;3 are based on edtimnated values of 
Cis, Crs, and Cys. 


J. W. SHELDON 
Computer Usage Company 
New York, New York 


B. ZONDEK 
Naval Proving Ground 
Dahlgren, Virginia 


M. FRIEDMAN 
Institute of Flight Structures 
Columbia University 
New York, New York 
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Optimum Quadrature Formulas In s Dimensions 


Although numerical procedures for functions of more than one variable are 
of considerable practical importance, they have received relatively little attention. 
So far as numerical integration is concerned, aside from successive application of 
appropriate one-dimensional results, as discussed in standard texts, and a few 
isolated special results for two dimensions, there has been little systematic study 
of the problem. Tyler [1] has collected a set of formulas of up to seventh degree 
for integration over rectangular regions in the plane, and up to fifth degree in 
three dimensions, and also gives a (2s + 1)-point third-degree formula for s 
dimensions. Hammer, Marlowe, Stroud, and Wymore [2, 3, 4] have pointed out 
useful methods of extending available results for spaces of a given dimension to 
higher dimensions, and to regions which are related to the given region by affine 
transformations. They have also developed a number of formulas for simplexes, 
cones, and spheres in s dimensions. 

It is the purpose of the present paper to describe a general approach to the 
problem which may often yield useful results and to use it to derive formulas in 
s dimensions using (s + 1) points for second-degree accuracy, and 2s points for 
third-degree accuracy. This method is applicable to any region, although its 
advantages are most pronounced for regions with a center of symmetry, and in 
our detailed calculations we will limit ourselves to a hypercube with center at the 
origin. We will be concerned with formulas of the form 


f- ‘ f o(x™, see, x dx i a dx™ = > c:o(x™, see, x) _ Re 
9 i=1 
where R¢y denotes the error of the formula for a particular function g. A formula 
will be said to be of degree n if Rg is zero whenever ¢ is a polynomial of degree n 
or less in the s variables x. 
If a formula is of degree n, the c; and x; must satisfy the set of equations 


(2) = G Il (x) = f we II (x')dx' i —_— 


j=l 
s 
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for all sets of ~; with sum equal to or less than m. An nth degree formula in s 
dimensions involves (n + s)!/n!s! such nonlinear inhomogeneous algebraic equa- 
tions, and each solution to such a set leads to an acceptable mth degree formula. 
Unfortunately, the task of solving systems of nonlinear algebraic equations is a 
formidable one, and general solutions have only been obtained in the simplest 
cases, although special methods have led to single solutions for more complex 
cases [1]. 

An instructive formulation of the problem can be obtained by a change in 
viewpoint. Let us define a set of m X m diagonal matrices by 


(3) G = [Veda] 
(4) XQ = [xi Sq, ]. 


In terms of these matrices, (2) becomes 


(5) tr {G TT (K)"G} = Iq, os, ap 


i=l 


Although the solution of a set of matrix equations such as (5) is generally no 
easier than the solution of a set of nonlinear algebraic equations such as (2), this 
formulation does emphasize certain considerations which are generally used in- 
tuitively in solving (2). In the first place, the fact that the validity of a formula 
of the form of equation (1) is independent of the numbering of the points is 
reflected in the fact that the trace of a diagonal matrix is unaffected by inter- 
change of diagonal elements. If the region S is symmetric with respect to change 
of sign of one of the variables, this fact must be reflected by an invariance of (5) 
to a change of sign of X“. If the region S is symmetric with respect to permuta- 
tion of certain variables, J,,,, ....,, must be invariant to permutation of the cor- 
responding 1;. 

Let us now specialize this to the case of hypercubes with edge length two, 
and center at the origin. In this case, the J,,, ..., », are independent of the order 
of the subscripts, and 


0 if at least one n; is odd 
(6) Be ng Sop ‘nel ed. 


II (n; +1) . 


j=l 


For a second degree formula, therefore, we have the equations 


(7) tr {GG} = 2° 
(8) tr {GXG} =0 

’ 2 
(9) tr {GXOXG} = — by 
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These equations in the traces may be converted to vector equations if we intro- 
duce the m-dimensional column vector with all elements unity, ¢, and its trans- 
pose, «7. Then if we define the vectors £, £;, &, ---, by § = Ge, &; = K%Ge, 
tj, = X®XGe and so on, (7), (8) and (9) become 


(10) e7GGe = (Ge)™(Ge) = E7E = 2° 
(11) e7GKGe = (Ge)?(K%Ge) = £7E; = 0 
(12) «™GXOX®Ge = (K%Ge)"(KGe) = £78 
= (K9XGe)"(Ge) = Eat = = bn. 


These, however, will be recognized as merely orthogonality relations aniong 
(s+ 1) vectors, &, :, ---,& and normalization requirements that |£|* = 2, 
|é;|* = 2*/3. Now (s + 1) orthogonal vectors span a vector space of dimension 
(s + 1) and this space must be a subspace of the vector space of dimension m 
consisting of all m-dimensional vectors. Thus, a second degree formula of the 
form of (1) can be obtained with m = s + 1, and for any higher value of m. 

Our argument has also furnished an explicit algorithm for constructing ex- 
amples of such formulas by orthogonalizing any linearly independent set of 
(s + 1)(s + 1)-dimensional vectors, and applying the proper normalization con- 
ditions. For example, if we orthogonalize the set, (1,1,1), (3, — V3 tan dv, 
V3 tan 8) and (3V3 tan #, 0, 0), in that order, we find 


(13) t = (2/V3, 2/V3, 2/3) 
2v2 2v2 Qe\ 2V2 4x 
(14) bs = (=> cos 9, “cos (0 + F), cos (o + F)) 
b= ( 3 ) sin 3 sin ( ;): 3 sin( 7)) 


so that our integration formula becomes 
1 1 
(16) f f ¢(s, t)dsdt 
-l¢’-1 
4 ’ 2x ae 2x 
- 3 { #(%F cos 2, V2 sin dv) + (VE cos(a + 7) V2 sin (2 _ +)) 


+ ¢ (Vacos(o +F), sin (a +F))| — Reg 


where Reg vanishes if g is a polynomial of degree 2 or less. The parameter # is 
arbitrary, although because of the periodicity of the trigonometric functions it 


may be taken to lie in the range — : <8< ‘ The vertices of any equilateral 
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triangle inscribed in a circle of radius v3 are thus seen to afford a satisfactory set 
of second-degree integration points. 

For higher degree formulas, the orthogonality conditions no longer suffice to 
define the minimai solution completely, but they still afford a substantial simpli- 
fication, as can be seen from the following discussion of third-degree formulas. 
For these, we have, along with (10), (11) and (12), the condition 


(17) e?GXOXMXGe = (KOXG6)7 (KG) = Eat = 0. 


Thus, we must consider the s(s + 1)/2 new vectors £, in addition to £ and the &;. 
These new vectors fall into two classes: the s £;; are orthogonal to every £:, but 
not to ~, while the s(s — 1)/2 &(7 # k) are orthogonal to both sets, or else are 
null vectors. 

Let us consider first the case where all the £ are null vectors. This implies 
that unless one or more elements of £ is zero, in which case the basic integration 
formula includes redundant points with zero weight, only one of the X“ can have 
any given element different from zero. The £;; cannot be null vectors, in view of 
(12), while from (11), unless the c; differ in sign, X“ must include elements of 
both signs, and thus at least two non-zero elements. We therefore conclude that 
the dimension of X“ must be at least 2s. 

For an equally-weighted formula, G is a scalar matrix, which we may write 
as gE, where E is the identity. From (7), for a 2s-point formula, g? must have the 
value 2*-'/s. For the minimum number of points, X will have but two non-zero 
elements of opposite sign, and, from (8), of equal magnitude, and these may be 
arranged in order such that 


(18) (X),; = x (Gni[s, 5-1 — 44, 5; ]) 


a diagonal matrix with the (27 — 1)th element equal to x and the 2jth to —x. 
From (9), it follows that 


(19) 2x'*g? = 2*/3 
so that for all j, 
(20) x@ = Vs5/3 


and we have the family of 2s-point third-degree integration formulas 





1 1 Qe-1 2s 
(21) f ee f o(x, sey x) dx oes dx) = yo o(x™, tee, x;) + R3¢@ 
vm | -1 i—1 


s 
where 


(22) xP = NE (5;, 25-1 — 54, 25). 


The case for s equal to three is included in Tyler’s collection [1]. For s greater 
than three, these formulas have the drawback that they depend upon values of 
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the function outside the range of integration, and thus may have larger remainder 
terms than are desirable. 

If one or more of the £ is non-null, there must be s + 2 or more orthogonal 
vectors, and the minimum value of m must be at least s + 2. This lower bound 
can be attained for s equal to two. In this case, we have a basis of four orthogonal 
vectors, £, £1, £2, and £12. If a four-point formula exists, these vectors completely 
span the space, and we can express £1; and £22 as linear combinations of £ and £12, 


(23) £55 = af + dk 


since £1; and £22 are orthogonal to £; and £2. If we multiply both sides of (23) by £7, 
it follows from (12) that a; must be 4. Hence, introducing the definition of £12, 
and rearranging, 


(24) XO(KO — oXM)¢ = be. 


Thus ~ is an eigenvector, and 4 an eigenvalue of the two matrices, 
X®(X™ — 5,X®) and X®(X® — b.X™). Since these matrices are diagonal, all 
the elements corresponding to non-zero elements of — must be equal to }. In 
particular, for equally-weighted four-point formulas, the diagonal elements of 
X® and X® must be the roots of the two equations 


(25) x)? — byxYx® = 8, x2)? — box Dx® = 3, 
which become 
(26) xO = (x — $)/dyx, 
and 

2 — bib2 + 5? 1 
27 ee Me, GE 
(27) ‘ 3 —bibs) * + 9 — bids) 


Now the sum of the squares of the four roots of (27) is equal to —2 times the 
second term, and if these roots are to satisfy (9) 


(28) 2(2 = bibs + b;*) - 3(1 _ bib2)4 
so that 
(29) by =. be = b. 
Accordingly, 
1 
30) 04 cs Bit? bs ee ee 
(30) Se +9 + 


(31) x = + {t(1+.-*_)} 
3 1+8 
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(32) xO = 4 {*(1 =*.)}' 


Thus, again we find that the integration points lie on a circle of radius V2/ a, 
and that there is a whole family of four-point third-degree cubature formulas, 
corresponding to the corners of squares inscribed in this circle. In three dimen- 
sions, it can be shown that a five-point third-degree formula is impossible, but 
that the vertices of any regular octahedron inscribed in a sphere of unit radius are 
satisfactory points for an equally-weighted quadrature formula. Even with matrix 
notation, the calculations for spaces of higher dimension, or for formulas of higher 
degree become extremely tedious. However, it seems likely that, as was observed 
for second- and third-degree formulas, the minimum number of points necessary 
for a quadrature formula of a given dimension, s, will continue to be of the order 
of some power of s, and not of the sth power of the number of points for a one- 
dimensional formula. 

It should be noticed that the application of the matrix-vector formulation is 
not limited to the hypercubical regions with constant weighting function which 
have been discussed in detail. It is equally helpful in deriving multidimensional 
analogues of the Gauss-Laguerre and Gauss-Hermite one-dimensional formulas, 
although the corresponding relations among the vectors are no longer simple 
orthogonality conditions, but require less obvious methods of solution. 

Because of the variety of formulas discussed, the problem of the truncation 
error for formulas of various degrees and dimension cannot be considered in 
detail. The more familiar methods for obtaining remainders in one dimension 
cannot, unfortunately, be extended to multidimensional problems, and Sard’s 
extension [5] of Peano’s theorem is one of the very few useful bases for an error 
estimate. Even so, the expressions are complicated, and depend upon several 
partial derivatives, thus requiring quite detailed information about the behavior 
of the integrand, which may often be unavailable. 

It is a pleasure to acknowledge the advice, suggestions, and encouragement of 
Dr. Gertrude Blanch, Aeronautical Research Laboratory, and Professor Thomas 
Holyoke, Miami University, Oxford, Ohio. 


Henry C. THACHER, JR. 


Aeronautical Research Laboratory 
Wright Air Development Center 
Wright-Patterson Air Force Base, Ohio . 
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TECHNICAL NOTES AND SHORT PAPERS 
A Method to Investigate Primality 


The method determines the smallest odd prime factor of a number N by 
testing the remainders left after division by the successive odd numbers 3, 5, - - - 
fm — 2, fm: here, fm is the largest odd number not exceeding N. If none of these 
remainders vanishes, N is a prime number. 

Let f be one of the odd trial divisors. Remainder ro and quotient go are defined 
by the relations 


N=10+ fao, 0< 10 < f. 
Now Qo is divided by f + 2, giving 
go = 11 + (f + 2)qi, 0<n<f+z2. 


Then q: is divided by (f + 4), etc., and this process is continued till a quotient 
(gn, SAY) equal to zero is found; r, is the last remainder in the sequence unequal 
to zero. After elimination of the g; we get the relations 


(1) N=rot+fnt+f(f+2)re+f(f+2f+4)rst+ --- 
+ f(f +2) --- Uf + 2n — 2)r, 
and 


(2) O<r <ft+2i. 


Once the sequence r; is known for a given value of f, it is easy to compute the 
corresponding sequence r;*, defined by the relations (1) and (2) with respect to 
f* = f + 2, as they are expressed in terms of the 7; by the recurrence relations 


(3) bo = 0, re = 7; — 206 + I)rigs — B+ Cf* + 22)din1, (¢ = 0,1, ---, m). 


The relation corresponding to (1) is satisfied for arbitrary values of the num- 
bers 6; with 7 > 1; they are fixed, however, by the relations corresponding to (2) 


(2*) O<r* < f* + 21. 

On account of the inequalities (2) and (2*)—and bo = 0—the 5; satisfy the 
inequalities 

(4) 0 <b; < 21. 

We have chosen by = 0. Then the relations (3) and (2*) with i = 0 determine 
ro* and bi; once 5; is known, (3) and (2*) with 7 = 1 determine r;* and dz, etc. 
The process is easily programmed. 


As fai = 0, and the inequalities (2*) with i = m are always satisfied with 
bn: = 0, the process terminates with 


r,* = 7, — ba. 
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As soon as 7,* = 0 is found—in that case it can be proved that r*,_1 ~ 0—the 
index m, marking the last r; ~ 0 in the sequence, is lowered by 1. 

In order to find the smallest odd prime factor of N, the 7; defined by (2) and 
(3) and f = 3 are computed. Here the only divisions in the process are carried 
out. At the same time the initial value of m is found. If N is large, this value may 
be considerable: for instance m = 11 is found for N ~ 10". The amount of work 
involved in each step is roughly proportional to m*. Fortunately large initial 
values of m decrease very rapidly. As soon as f-(f + 2)-(f +4) > N, m takes 
the value 2. This is its minimum value: when 7,* = 0 with = 2 is found, 
(f* + 2)? > N holds and N is a prime number. (If not, we should have found 
an 79 = 0 earlier and should have stopped there.) 

The process still may be speeded up. Let 0d,’ be the minimum of 5, for fixed 
n up till a certain moment: then it can be shown that the next 6, satisfies 


(5) b, <b,’ + 1. 


Let us apply this to the last stage m = 2. According to (4) bz satisfies 0 < be < 4. 
According to (5), however, the only possible values for bz are 0 and 1 as soon as 
a value b: = 0 once has been found. This is bound to happen for f ranging 
(roughly) from (4N)? to (8N)%. In the case b2 = 0 it is apparently unnecessary 
to test whether 72 = 0 is reached. (If N > 144, the case b, = 0 with m = 2 occurs, 
before 7,* = 0 with m = 2 is found; prime numbers are then always detected in 
this last stage.) 

The less efficient steps of the process for large m (i.e., small f) could be avoided 
by carrying out divisions for small values of f (see Alway [1]). However we 
strongly advise against doing this. . 

If the process described above is started at f = 3, the whole computation can 
be checked at the end by inserting the final values of f and 7; into (1). As all the 
intermediate results are used in the computation, this check seems satisfactory. 

If a double-length number N is to be investigated, another argument can be 
added: division of N by small f may give a double-length quotient, i.e., two 
divisions (and two multiplications to check) are needed for each f. In our case 
only part of the initial m divisions are double-length divisions. 

The process described above has been programmed for the ARMAC (Auto- 
matische Rekenmachine van het Mathematisch Centrum). The speed of this 
machine is about 2400 operations per second. A twelve decimal number was 
identified as the square of a prime in less than 23 minutes. 


E. W. DijyKsTRA 


Mathematical Centre 
Computation Department 
Amsterdam, The Netherlands 


i. G. G. Atway, “A method of factorisation using a high-speed computer,” MTAC, v. 6 
1952, p. 59-60. 
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Equally-Weighted Quadrature Formulas 
for Inversion Integrals 


In a previous article [1] the author considered Gaussian-type quadrature 

. ‘ eae We es 
formulas for the numerical evaluation of inversion integrals ar - F(p)dp, 
mt 


where an ”-point formula was exact whenever F(p) was a polynomial of degree 
(2n — 1) in 1/p. In this present note we consider equally weighted (i.e., Chebyshev 
type) quadrature formulas of the form 


1 io ep a 
(t) Pride un p POMP = 5% Fibs) 


where (1) is exact whenever F(p) is any polynomial of degree m in 1/p. The 
analogous set of equally weighted quadrature formulas for evaluating infinite 
integrals that are direct Laplace transforms has already been considered in one 
of the author’s earlier papers [2]. Similar to the derivation given there, it is 
easily seen here that the well-known relation 


-. ety a _ 4 
Qrid-in p\p ee. 


by choosing r = 0, establishes the factor 1/m outside the summation in (1), and 
the choice of r = 1, 2, ---, establishes the following necessary and sufficient 
conditions on ;, in order that (1) should hold whenever F(p) is an arbitrary nth 
degree polynomial in 1/p 


(2) (2) =5, poor © Tee 


j= 


Following the usual methods [2], one determines for the m-point case of (1) the 
coefficients of the polynomials ¢,(z) whose zeros 2; are the reciprocals of the 
required points p;. Thus from ¢$,(z) = 2* + a2"! + aoa"? + --- + Gp_id + Gn, 
the coefficients a, are found successively from 


(3) Ray + Gy-1 dX 2; + axe 27 + +++ $ai Def" + D 2# = 0, 


i=l j=l i=l i=l 


k = 1,2, ---,m. 


Below are the general formulas for the first ten coefficients @;, @2, - - -, @i9 for any , 
having meaning, of course, only for an a, where k < n 
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aq, = 


ae 


ag 


ay = 


a10 


n* n* 25 


a— ap an one of 


24. +8 288” 














n 



































6 ’ 
ba eee Bip sy ies % 
~ ~ 720" 24°” 288” * 288 600’ 
n® n® 43 25 1507 n 
Se sn —_ ate a ae ae 
720 96+ 1728” ~ 1152” * 259200" ~ 4320’ 
Te 13 1741 53 s 
= es ‘teem 5 S at 3 eee ee ® 
5040 480 1728” * 1152” 250200” * 43200” ~ 35280’ 
n’ n? 61 7 17627 3779 
= ae © tak deci 5 yy 3 
40320 2880 * 34500” 1728” * 4147200” — 20 73600” 
15787 : n 
n —_— 
677 37600 3 22560’ 
n? n’ 7 19 22289 
-_, lie 7 cco <eneae 5 
3 62880 * 20160 20736” * 17280” ~ 124 41600” 
ea 12317, " 
— 98 — eee 99 n? — ; 
20 73600 10973 49120 3048 19200 32 65920 
ni a ae Sa 
= - n> — n™ + ——_____ 
36 28800 1 61280 14 51520 46080 1244 16000 
11347, 22 37339, 26791 


~ 165 88800” — 54867 45600” 2709 50400” 


5.90383, - 
9 14457 60000” — 362 88000 





The first ten polynomials ¢,(z) are 


¢1(z) =.= a 


$2(z) = 2? — +5, 


15 
(2) = 2 — 32° +7 2-75 


29 
12’ 





Oo Sem Fe aes oc oSlhUutlCUCO = eal 


oc @®D 






20’ 
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g(a) = st — 48 + 78 — SS, 

os(a) = a6 — Set Se — Pe 4 EO 

62) = - 68 +S — Se 4 Se ST, ae 
or(a) = a? — Tat pat — gy SP Bi 


53 02619 180 44381 


1 29600 ° 9 07200 ’ 


628 4037 277 1 
os(z) = 5 — 857 + 30s* — —— a6 + 3 3 9 22919 





36 25 8100 


29 22187 . 1455 11171 
39690 42 33600 ’ 


153 407 3027, 1 03911 6 50239 


























ia Of bee in co 
o(s) = 2 — Sst +s 1* ts"? 400 +00 7 
_ 84 99571 414 78457 15146 11753 
39200 3 13600 254 01600 ’ 
95 1280 . 43235 1 70381 14 90251 
— Oe Se Pe ii 5 
tel a - SAR EA 9 * +14" 300.” 2502” 
173 63761 4151 $8089 , 32552 25203, 14 23249 22009 
31752 10 16064 137 16864 13716 86400 


The values of »; and the zeros 2; = 1/p; of the above polynomials ¢,(z), for 
n = 1(1)10, are given in Table 1. 

The zeros 2; = 1/p; of ¢,(z) were calculated for m = 3,4, ---, 10, by first 
obtaining an initial approximation, using a procedure that had been employed 
upon the Univac Scientific Computer (ERA 1103) at the Convair Digital Com- 
puting Laboratory. The initial approximations to the complex zeros were then 
used to construct approximate real quadratic factors, which were refined by 
Bairstow’s method (Milne [3], Olver [4]), using only a desk calculator. The 
initial approximations to the real zeros were refined by Newton’s method. All 
factors of ¢,(z) were checked by four different formulas (see [4], p. 414). Also a 
final functional check was performed upon the values of 1/p; by substituting into 
equation (2) for m = 1(1)10, and r = 1(1)m. 

The 8-decimal values of »; and 1/p; in Table 1 are guaranteed as far as the 
seventh decimal place. But they are believed to be correct to within around two 
units in the eighth decimal place for m = 1(1)7 and have a high probability of 
being correct to within several units in the eighth decimal place for m = 8(1)10. 
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TABLE 1. p; and 1/p; 


n j D; 1/p; 
1 1 1.00000 000 + .00000 000z 1.00000 000 + .00000 000: 
2 4.2 .66666 667 + .47140 4527 1.00000 000 + .70710 678 
3 2 .46343 318 + .66891 6552 -§69981 792 + 1.01011 279 
3 .62485 778 + .00000 0002 1.60036 417 + .00000 000: 
4> 4,3 .31209 699 + .78442 8702 .43788 772 = 1.10059 2777 
3,4 .54603 449+ .22670 4977 1.56211 228 + .64856 456: 
5 PO .19029 304+ .86260 4997 .24387 201 + 1.10548 034: 
3,4 46724 697 + .36843 4487 1.31966 923 + 1.04058 814: 
5 .53392 634 + .00000 000z 1.87291 752 + .00000 000: 
6° 1,2 .08786 626+ .92009 404: 10285 254 + 1.07702 3312 
3,4 .39416 727 + .46819 7997 1.05229 916 + 1.24993 725: 
5,6 .49826 825+ .14769 9207 1.84484 830 + .54685 928: 
7 1,2 —.00076 496 + .96470 8252 —.00082 196 + 1.03658 217: 
3,4 .32727 973 + .54346 944: -81317 581 + 1.35033 175: 
5,6 45588 935 + .25464 1182 1.67190 107 + .93385 56% 
7 49224 949 + .00000 0002 2.03149 016 + .00000 000: 
8 2 —.07902 919 + 1.00066 4802 —.07843 500 + .99314 110: 
3,4 .26601 917 + .60293 7622 61252 402 + 1.38829 762: 
5,6 41223 2514+ .33698 9857 1.45409 421 + 1.18868 5952 
7,8 47182 912 + .10911 5332 2.01181 677 + .46525 329% 
9 1,2 —.14919 526 + 1.03046 7527 —.13761 844 + .95050 835: 
3,4 .20966 304+ .65149 3532 44761 307 + 1.39088 424: 
5, 6 .36931 455+ .40305 3927 1.23580 348 + 1.34870 244: 
7,8 44525 659+ .19444 9152 1.88616 975 # .82371 403: 
9 46815 071 + .00000 000 2.13606 428 + .00000 000: 
10 M —.21284 773 + 1.05570 9532 —.18351 682 + .91023 036 
‘ 15754 418 + .69213 4697 .31266 794 + 1.37363 580: 


32790 360 + .45764 0257 1.03454 187 + 1.44386 337: 
41610 417 + .26374 9507 1.71443 371 + 1.08670 1517 
45488 509+ .08636 297: 2.12187 330 = .40285 183i 


o 


- 


Oo NTN 
— CO CO ® NO 


The refinement of the factors and zeros of ¢,(z) by Bairstow’s method was 
done by Mrs. Genevieve Mullin Kimbro. 


HERBERT E. SALZER 
CONVAIR Astronautics 
San Diego, California 


1. H. E. Sauzer, “Orthogonal polynomials arising in the numerical evaluation of inverse 
Laplace transforms,” MTAC, v. 9, 1955, p. 164-177. 

2. H. E. Satzer, “Equally weighted quadrature formulas over semi-infinite and infinite 
intervals,” Jn. Math. and Physics, v. 34, 1955, p. 54-63 

3. W. E. MILne, Numerical Calculus, Princeton University Press, Princeton, New Jersey, 
1949, p. 53-57. 

4. f. W. J. Otver, “The evaluation of zeros of high-degree polynomials,” Roy. Soc. London, 
Phil. Trans., v. 244, Ser. A, 1952, p. 385-415. 
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An Open Formula for the Numerical Integration 
of First Order Differential Equations 


1. Introduction. The integration scheme given below combines the accuracy 
of multi-point formulas with the convenience, lack of starter formulas, etc., of 
two-point formulas, as described, for example, by Milne [1]. The price paid is 
the necessity of evaluating the right hand side of 


(1) y’ = f(x,y) 


at points which may lie outside the range of integration. It is the writer's belief 
that this is compensated for by the ease of programming the scheme for a digital 
computer. 


2. The Integration Formulas. For each n = 3,4, 5, ---, there exists a matrix 
D,, of order m and rank m — 1 such that if 9 = (yo, 1, ---, ¥a-1) is a vector of 
values of the function y(x) at equally spaced points xo, x9 + h, ---,xo + (m — 1)h, 
then 9’ = (yo, v1’, «+ -, ¥’n—1) is given by 
(2) 7 =Dg+E 


where E£ is an error vector whose elements are O(h*-1y™ (é)). 
Recurrence formulas for the propagation of the solution of (1) can be obtained 
by setting 


(3) Dg + E = f(z, 9) 
where 


f(Z, 9) = (f (xo, yo), f (x1, y1), a f (%n-1, Ya—1)); 


from which all but two variables can be eliminated successively, yielding the 
desired recurrence relation. Thus for » = 3, 


1 h? ; 
(4.1) yo’ = — (— 3y0 + 41 — 92) += y™ = fF (Xo, Vo) 
2h 3 
Lire h? 
(4.2) vy = ah (— yo + ¥2) — ae = f(x1, y1) 
aot he 
(4.3) yy = (yo — 4y1 + 3y2) + ; af = f (x2, y2). 


Solving (4.2) for ye, inserting in (4.1) we get 
h h* 
(5) v1 — Yo= 2 {f (xo, yo) + f(x1, 91)} — 2?” 


the familiar trapezoidal rule. 
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With n = 4 we eliminate ye, y; from 


1 8 
(6.1) yo = oh (— 11y0 + 18y1 — 9y2 + 2ys) — li = f (xo, Yo) 
i aong y 
(6.2) no (— 2¥0 — 391 + 6y2 — ys) + 73% > = f (x1, y1) 
—— ae 
(6.3) rs) (yo — 6y1 + 3y2 + 293) — nD” > = f (x2, 2) 


1 8 
(6.4) = oh (— 2y0 + Dy: — 18y2 + 11ys) + 2" = f (Xs, ys) 
getting 


h 4 
(7-1) ¥1 — Yo = Fy [SF (%o, yo) + BF (%1, 91) — F (x2, ¥2)] + a 


h* 
(7.2) y2 = Syo — Ay + 2h{ f (xo, yo) + 2f (x1, ¥1)} + “7 
Replacing y2 in (7.1) by (7.2) and using the mean value theorem we get finally 
h 
(8.1) 1 — Yo = Fy {SF (%o, Yo) + Bf (x1, 1) — f(x2, ¥2*)} 


O - y¥9 (2) [1 +o ee, 0| 


3 dy 
where 
(8.2) yo* = Syo — 4y1 + 2hLf (xo, Yo) + 2f (x1, 91) ] 
(8.3) xoSES x 


and 7 lies between ye and y2*. 

If f(x, y) is linear in y then y is given explicitly by (8.1) and (8.2) ; otherwise 
an initial guess may be made of ¥;, y2* computed from (8.2), and a new ¥:, com- 
puted from (8.1), repeating the process until convergence is reached before pro- 
ceeding to the next point. 

In the linear case, (1) takes the form 


(9) y’ (x) = P(x) + Q(x)y(x), 
and (8.1), (8.2) become 


‘ 


10) [1-4 2a+09 + 0.0] -»[1+4 @- 0 -*00,| 


h h? h* h 
<< [S5Po + 8P; — P2] — = [Po + 2P1]+ m2 @) [1 +a]. 
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The method clearly generalizes to n > 4 yielding formulas of higher order 
accuracy. The involved calculations necessary in carrying out the propagation, 
however, may make the method too inconvenient for machine calculation. 

3. An example. We illustrate the method with the solution of 


F mids 
(tt) pi 2 


for which the exact solution is y = 3e* — 1. 
The functions P(x), Q(x) are here identically unity, and (10) becomes 


h2 h? h? 
(12) [1:-a+%]y,=-[1-*],4+0-% 
or 
(13) yi = A(1 + yo) — 1 
where 


= (1 - wl — a+ ep. 
The fourth order Runge-Kutta method also gives (13) with 
= (1 +h+ 3h + 3h’). 
For h = .05 the solution was propagated twenty steps with the following results: 


x y (Runge-Kutta) 


y (Eq. 10) y (Exact) 
0.00 2.00000 0000 2.00000 0000 2.00000 0000 
0.20 2.66420 4603 2.66420 4351 2.66420 8275 
0.40 3.47546 5123 3.47546 4508 3.47547 4093 
0.60 4.46633 9967 4.46633 8842 4.46635 6401 
0.80 5.67659 6021 5.67659 4190 5.67662 2786 
1.00 7.15480 4623 7.15480 1828 7.15484 5486 





4. Significance. In order to obtain accuracy of higher order than the third, 
using the Runge-Kutta technique, it is necessary to evaluate the functions Q(x), 
P(x) in equation (9) at points interior to each mesh interval. On the other hand, 
equation (10) requires that evaluation only at points at which y is to be calculated, 
plus one extra beyond the right hand limit of integration. 

This fact, plus the freedom from starter formulas required by some multi-point 
methods combine to enhance the value of equation (10), or equivalently (8.1), 
(8.2) for machine computation. 


HERBERT S. WILF 


Nuclear Development Corporation of America 
White Plains, New York 


1. W. E. Mitne, Numerical Calculus. Approximations, Interpolation, Finite Differences, 
Numerical Integration, and Curve Fitting, Princeton Univ. Press, New Jersey, 1949, p. 96-98. 
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Note on “A Method for Computing Certain 
Inverse Functions” 


The method for computing certain inverse functions, one binary digit at a 
time, which was described by D. R. Morrison in a recent issue of MTAC [1] 
has been used in this laboratory. In particular, a subroutine for computing an 
inverse cosine, x = arc cos y, based on the method was given in [2] (part III, 
subroutine T4, p. 152-153). It was, however, pointed out by van Wijngaarden [3] 
that the method gives poor accuracy for certain values of the argument, namely, 
those for which one or more of the functions cos x, cos 2x, cos 4x --- cos 2*x are 
near unity. When x is near zero the error is, perhaps, of little importance since the 
equation x = arc cos y does not then determine x with any great precision, but 
this is not the case when «x is near 1/2, r/4, etc. In general abnormally large errors 
may occur if, in Morrison’s notation, 


dy,/dx = 0 for any n < N, 


since 6 will then be of order ve if d?y,/dx? + 0, and of larger order otherwise. The 
number of correct figures in the value obtained for an inverse cosine or similar 
function may, as a result, be only about half as many as Morrison suggests. 

Although, in cases in which the above objection does not apply, digit by digit 
methods of computing functions may sometimes be useful in a digital computer 
on account of their simplicity, they are, in general, slow in operation, and unless 
storage capacity is very restricted other methods are, in our experience, generally 
to be preferred. 


M. V. WILKES 


D. J. WHEELER 
University Mathematical Laboratory 


Cambridge, England 


1. D. R. Morrison, ‘‘A method for computing certain inverse functions,” MTAC, v. 10, 1956, 
p. 202-208. 

2. M. V. WiLxkgs, D. J. WHEELER, & S. GILL, The Preparation_of Programs for an Electronic 
Digital Computer, Addison-Wesley Press, Cambridge, Mass., 1951. 

3. A. VAN WIJNGAARDEN, “Erreurs D’Arrondiment dans les Calculs Systématiques,”’ Centre 
National de la Recherche Scientifique, Colloques Internationaux, v. 37, 1953, p. 285-293. 


REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


75[A ].—WALTER Scumipt, Der Rechner, Technischer Verlag Herbert Cram, 
Berlin, 1955, xix + 200 p. DM 18.00. 


This gives on page m, m = 1(1)200, the product mn, to one place of decimals, 
where m = p + 6, p = 0(1)100, 6 = 0.1(.1).9; 2, 4, 2, 2, 4, 2. The format is 
satisfactory and the printing tolerable. The accuracy has not been checked. 

There is a rather lighthearted introduction, which gives various examples of 
the use of the table. 
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This review was prepared by J. Todd for Mathematical Reviews. 


76[A ].—FRANZ TRIEBEL, Rechen-Resultate. Achte Auflage, Technischer Verlag 
Herbert Cram, Berlin, 1956, ii + 285 p. DM 26.00. 


This gives, on page m + 2 for m = 1(1)100, the product mn to two decimal 
places, where m = 1(1/4)100. There follows the product mn for m = 101(1)1000, 
n = 1(1)100, each page covering five values of m; there is also a table of mn, 
m = 1(1)300, » = 3, 3, . 

There is a brief introduction showing the use of the tables. There are elaborate 
thumb indices. Printing is satisfactory; the accuracy has not been checked. 

eS aA 


This review was prepared by J. Todd for Mathematical Reviews. 


77(B ].—H. NacGuer, Table of Square Roots of Integers, on microfilm, 101 frames, 
deposited in UMT FILe. 


This table gives VN to 15D (17S) for N = 1(1)10,000. The table was com- 
puted on an Elliott 402 digital computer, and printed on an automatic printer 
from five-hole punched paper tape. The author states that the method of compu- 
tation insures that the error is within one unit in the last decimal place. A random 
spot check by the reviewer with 18D hand computed values revealed an error of 
1 in the last place in two cases out of ten, the other eight cases being correct. 


Davip A. Pope 
University of California 
Los Angeles, California ; 


78[B, C, D, E].—V. BELevitco & F. Storrer, “Le calcul numérique des fonc- 
tions élémentaires dans la machine mathématique IRSIA-FNRS,”’ Acad. r. 
de Belgique, Ci. d. Sciences, Bull., s. 5, v. 52, 1956, p. 543-578. 


This article treats in detail the problems of approximating elementary func- 
tions by polynomials, for use on a digital calculator using floating decimal arith- 
metic. The specific calculator the authors have reference to is the Belgian machine 
ISRIA-FNRS [1], which uses a floating decimal arithmetic with 15 significant 
digits in the mantissa, and an exponent between — 50 and 50. The basic command 
structure allows the machine to add, subtract, and multiply; all other arithmetic 
operations must be programmed. 

It is clear that in dealing with floating point operations, the desirable criterion 
for an approximating polynomial is that the relative error, rather than the absolute 
error, of the approximation be bounded. Two methods of deriving a polynomial 
with a given relative error from a polynomial with a given absolute error are 
discussed ; the first method depends on the property that if the function f(x)/x 
is approximated by a polynomial p(x) on an interval J with an absolute error = e, 
then the approximation xp(x) to the function f(x) will have a relative error = Me, 
where m = sup|x/f(x)|, for x in J. The second method utilizes the fact that if 
the derivative f’(x) is approximated by a polynomial with a certain absolute 
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error, then the integral of the polynomial will approximate f(x) with a bounded 
relative error. 

Several methods of obtaining the approximating polynomials are discussed, 
and illustrated by deriving the polynomials used in the calculation of 1/x, x7, 
sin x, arc tan x, 10°, 10* — 1, logio x. Tables are given of the coefficients of the 
polynomials derived, to 15S. 

An analysis of the rounding error in the calculation of polynomials is also 
made, and the question of which elementary functions should be chosen as basic 
is investigated. For example, for x near zero, the function 10” — 1 is chosen as the 
basic function, and 10? is gotten from it by addition of unity; this has obvious 
advantages in the calculation, for example, of sinh x for x near zero. 

The basic law of floating point coding is ‘‘avoid subtracting two nearly equal 
numbers,” and the authors have shown ingenuity in complying with this law 
while calculating the elementary functions. The ideas contained in this article 
should prove valuable to everyone who is coding in floating decimal (or floating 
binary). 

Davip A. Pope 


University of California 
Los Angeles, California 


1. OrFicE oF Nava RESEARCH, Department of the Navy, A Survey of Automatic Digital 
Computers, Washington, D. C., 1953, p. 53. 


79[D].—L. S. KureNov, Pyatiznaisnye tabliisy trigonometricheskikh funktsii s 
argumenton, vyrazennym v isasovoi mere (Five-place tables of the trigonometric 
functions with argument expressed in hourly measure). Second edition. Gosu- 
darstv. Izdat. Tehn.-Teor. Lit., Moscow, 1954, 172 p. Price 7.70 rubles. 


The main table gives values of the six trigonometrical functions at interval 
of 4 seconds of time (i.e., 1 minute of arc), usually to 5S, in the range 0 to 3 
hours. (Whenever the leading figure is unity, five further digits are given.) First 
differences are usually given, with their proportional parts alongside the tables. 
An auxiliary table, new in this second edition, gives cot x, cosecx to 5S for 

= 0(.1°)8™(1*)40™. There is a table of sin? (x) for x-= 0(4*)12". Entries with 
terminal 5 are marked with + to indicate in what direction rounding should be 
made. There is a collection of conversion tables, constants and formulae from 
plane and spherical trigonometry and there are worked examples showing the 
direct and inverse use of the tables. The tables are clearly printed. 

There are no references to sources, nor is there a description of the construc- 
tion. It is stated that the entries are correct to half a unit in the last place. The 
tables should be convenient for those who require something between the two 
volumes issued by L. J. Comrie [1] which gives 4+ decimals at 10* interval, and 
the British Nautical Almanac Office [2] which gives 7D at 1 interval. 

te 

This review was prepared by J. Todd for Mathematical Reviews. 

L. M. MiLne-THomson & L. J. Comrie, Standard Four-Figure Mathematical Tables, Edition 
B, MacMillan & Co., Ltd., London, 1931; Edition A, MacMillan & Co., Ltd., London, 1944. 


H. M. NavticaL ALMANAC OrrIce, Seven- - Figure Trigonometrical Tables for Every Second 
of Time H. M. Stationery Office, London, 1939. 
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s0[E, H, P].—J¢rceEn RyBNER, Nomogrammer over komplekse hyperbolske funk- 
tioner (Nomograms of Complex Hyperbolic Functions), Jul. Gjellerups Forlag, 
Copenhagen, 1955, 39 + 60 p. of illustrations, 30 cm. Price Dan. Kr. 44.00. 


This is a second edition of this useful book. (The first edition was reviewed by 
K. G. Van Wynen, RMT 526, MTAC, v. 3, 1948, p. 174-175.) The unhandy 
binding of the first edition has been improved, errors have been removed, notation 
has been changed (complex arguments are now A + jB instead of b + ja), 
nomograms have been added extending the range of the argument for hyperbolic 
tangent, and nomograms for interaction loss and interaction phase shift have 
been added. 

In addition to discussion and formulas, alignment charts are included for 
sinh (A + 2B) and cosh (A + 1B), 0 <A <4, for tanh (A +1B),0<A<S3 
and for various incidental functions as follows: 


x + zy = r(cos@ + isin 8); 
R(cosa +isina) = 1+7r(cosé+isin@), r <1; 





4, = in(**2), i < |Z| < 10; 

B, = arg(*47), 

A, = $l1n (1 — 2e*4 cos2B + e“4), OS ASI, OF BK 9°; 
B, = arctan5>"— A>0; 

f=. 10-* < L < 10° henry, 10-" < C < 10 farad 


K = NS 10° << L < 105, 10° <C <10. 


The tables and their units are chosen for convenience in electrical filter design 
and other similar problems. Accuracy is adequate for most such applications. 

The author notes the following errata: 

“On the four new charts for tgh (A + jB) = 7/@ covering the ranges 
A = 2,00-2,25; 2,25-2,50; 2,50-2,75; 2,75—3,00 nepers, the B,---B, scales are 
erroneously marked ay: - «a. 

In the list of contents, page 6, the formula 7 for the interaction loss should read : 


A, = }1n (1 — 2e~*4 cos 2B + e~*4). 
In the Danish preface, page 9, line 7 from the bottom, the word mono- 


grammerne should read nomogrammerne.” 
©. &: F. 
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81[E].—Homer S. Pow ey, Table of log cosh x. One typewritten sheet, 28 cm., 
deposited in the UMT FILE. 


This table lists log cosh x, 3D, for x = 10(.5)20(5)40. 
Because of the high value of x, the table is essentially x log e — log 2. 
ca TF. 


82[F, Z].—CarL-Er1K FROBERG, Hexadecimal Conversion Tables, C. W. K. 
Gleerup, Lund, Sweden, 1957, 26 p., 22 cm. Price 3 kr. 


A revised version of the conversion table so widely used around computers 
with hexadecimal input. (See RMT 1042, MTAC, v. 7, 1953, p. 21.) In this 
edition the binary point follows the first binary digit of fractions; the first digit 
is used as a sign. This is clearly important in use of the table; for some machines 
the numbers must be doubled and the sign inserted properly. 

Using A, B, C, D, E, and F for the hexadecimal digits ten through fifteen 
respectively the table lists: 


1. Decimal and hexadecimal integers over the decimal ranges 1 (1)1024(16)4096 
and 10*(10*)10*+!, k = 2(1)12; 

2. Hexadecimal equivalent of decimal fractions x, x = 10-*(10-*)10-**, 
k = 2(2)16; 

3. Conversion of 2-10*, to normalized hexadecimal numbers, = 1(2)9, 
k = —12(1)12 and conversion of 10* and 10~ for k = 13(1)25; 

4. Hexadecimal form of constants frequently met; 

5. Decimal form of hexadecimal fractions x = 16-*(16~-*)16-*+! (subject to 
the sign convention mentioned earlier) k = 1(1)10. 


For similar octal-decimal tables, see van Wijngaarden [1], Karst [2], and 
Causey [3]. None seems to have been widely distributed. 
The following erratum was communicated to the author by S. Arkéus: 
On page 23, for "log 10 = 6A4D3 C25E8 1209D 82 read 6A4D3 C25E6 
8DCS58 82. 
i. a ws 


1. A. VAN WIJNGAARDEN, “‘Decimal-binary conversion,’’ Report R-130 of the Computation 
Department, Mathematical Center at Amsterdam (see following review 83). 

2. EpGAR Karst, Tables for converting 4 digit decimal fractions to periodic octal fractions. [See 
Review 6, MTAC, v. 10, 1956, p. 37.] 


3. Rospert L. Causey, Decimal to octal and octal to decimal conversion tables. [See Review 65, 
MTAC, v. 10, 1956, p. 227.] 


83[F, Z].—A. vAN WIJNGAARDEN, Decimal-Binary Conversion and Deconversion, 
Report R-130 of the Computation Department of the Mathematical Centre, 
Amsterdam, 1951, 41 p., mimeographed, 33 cm. 


This useful table, hard to read, gives decimal equivalents of octal numbers 
0(10)303230 (all digits OCTAL!). In addition it gives decimal values of 2", 
n = 1(1)50, exact, decimal values of 2—-", » = 1(1)50, 20D, octal values of 10, 
n = 1(1)18, exact, and octal values of 10-", m = 1(1)18, to twenty octal digits. 

References to similar tables and their usage may be found in Review 82, above. 
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s4[F, X].—P. Davis & P. Rasrnowitz, “‘Abscissas and weights for Gaussian 
Quadratures of High Order,”” NBS Jn. of Research, v. 56, 1956. 


This paper contains 20D values of weights and abscissas for Gaussian quadra- 
ture rules with m = 2, 4, 8, 16, 20, 24, 32, 40, and 48 points. Corresponding values 
are available from the authors for the cases nm = 64, 80, and 96. 

The weights, a;,, and abscissas, x;,, enter the approximate formula 


f "flelde ~ E draf (ian). 


k=1 


These numbers, which have been adequately checked by the National Bureau 
of Standards, should prove useful in cases where: (a) the integrand f(x) can 
feasibly be computed for arbitrary x and (b) a high degree of accuracy (result is 
exact if f(x) is a polynomial of degree < 2m — 1) is needed. 


T. H. SouTHARD 
University of California 
Los Angeles, California 


85[I, X].—K. A. Karpov, Tabliisy Koéffisientov interpoliaisionnot Formuly 
Lagranzha (Tables of Lagrangian Interpolation Coefficients), Akad. Nauk 
SSSR, Moscow, 1954, 79 p., 26 cm. Price (including [3 ]) 61 rubles. 


These tables contain for four-point Lagrangian interpolation the coefficients 
A,(é),i = —1(1)2,¢ = —1(.001)2 and for five-point interpolation the coefficients 
A,(t),i = —2(1)2, # = —2(.001)2, 6D. It was issued as a supplement to [3]. 

Several entries were checked against the Mathematical Tables Project’s more 
extensive tables [1], and no discrepancies were found. The printing of the present 
volume is clear and easy to read, and the small size of the volume makes the tables 
much handier than [1 ] for the many times when four- or five-point interpolation 
with these increments suffices. 

The numbers listed seem to have been rounded individually rather than as a 
group, so that 2,A,(#) may differ from 1 in the last digit. This happens for 
t= 0.001 in the four-point coefficients, for example, where the entries are 
A_, = —0.00033 3, Ao = 0.99949 9, A, = 0.00100 0, and Az = —0.00016 7. 
For a slowly varying function tabulated to many places subtraction of a constant 
(usually common leading digits) may be required to reduce the absolute value of 
the tabulated entries in order not to introduce a rounding error; restitution of 
the subtracted portion must follow the interpolation. This process is not always 
convenient, and the reviewer would prefer the forced rounding used in [1] to 
assure that £;A;(#) = 1. It would seem to be reasonable to mark each digit which 
should be rounded up in a further truncation of the coefficients so that this unit 
sum could be maintained ; thus if only 3D values are needed the user would round 
up an overscored third digit and leave others as printed in the table. The addi- 
tional cost of printing might well be justified in ease of using the tables. 

An estimate of the continuing need for Lagrangian interpolation is made in 
the MTAC review of [1]. More extensive tables are available in [1] (and this is 
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noted in the present work), and tables for sexagesimal arguments are available 
in NBS AMS No. 35, [2]. 
a 


1. NYMTP, A. N. Lowan, technical director, Tables of Lagrangian Interpolation Coefficients, 
Columbia University Press, New York, 1944. [RMT 162, MTAC, v. 1, p. 314-315.] 

2. NBS Applied Mathematics Series, No. 35, Tables of Lagrangian Coefficients for Sexagesimal 
~~ So U. S. Govt. Printing Office, Washington, D. C., 1954. [Rev. 54, MTAC, v. 11, 
p. 108. “a od 

3. K. A. Karpov, Tablitsy Funktsit w(s) = e~®* fo e** dx v Kompleksnoit oblasti, Akad. Nauk 
SSSR, Moscow, 1954. 


86[K ].—BeELt ArrcraFtT CorporaTIon, Table of Circular Normal Probabilities, 
Operations Analysis Group, Dynamics Section, Report No. 02-949-106, 1956, 
iv + 305 p., 22 X 27 cm. (oblong). A limited number of copies are available 
by writing to the Research Division, Bell Aircraft Corp., Buffalo 5, New York. 
One copy deposited in the UMT FILe. 


Circular normal probability integrals 


1 _ (e—a)?+(y—b)? 
P(T, «) = = - dyd. 
(T, o) oxo? a Yax 


for R/a = 0(.01)4.59 and Va? + &/¢ = D/o = 0(.01)3, SD. 

These tables were calculated on an IBM model 650 computer and checked 
against other available tables with perfect agreement over regions of overlap. 
Two thousand values chosen randomly were checked and first and second differ- 
ences were examined. 

The calculation is described and a most elementary illustrative example given. 

Printing by a photographic offset process is adequate. 

No differences or aids to interpolation are given. 

C. B. T. 


87[L].—T. Pearcey, Table of the Fresnel Integral to Six Decimal Places, Cam- 
bridge, England, at the University Press, 1957, 63 p., 24 cm. Price $2.50. 


This table was compiled for, and printed by, the Commonwealth Scientific 


and Industrial Research Organization, which published an Australian edition at 
Melbourne in 1956. The functions tabulated are 


1 1 “cost. * 
c=5f J4()dt = f Sa 
2J @) V(2m) Jo Vt 

= sin t 


1 1 
S=5f Aa = —— Sy 


Tables with argument u, where x = }2u?, are also common, but the argument 
of the present table is x, though a few expansions in terms of u are displayed in 
the explanatory text. The values of C and S are listed to 7D for x = 0(.01)1 and 
to 6D for x = 1(.01)50, with 6 throughout. 

The values were mostly computed by subtabulation of a 7D table for 
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x = 0(.02)1 given by Watson [1 ], and a 6D table for x = 1(.5)50 due to Lommel, 
and reproduced by Watson ; Lommel’s table is unreliable in the sixth decimal, but 
was differenced and corrected by recomputation where necessary. The compiler 
states that “errors in rare cases may amount to 1 unit in the last figure.” 

The reviewer compared the values up to x = 2 with 8D values (with possible 
error up to 2 final units) computed by Corrington [2]. Apart from a number of 
cases (almost all in the 7D portion) in which the two tables differ by between 
0.5 and 0.9 units of the last place printed in the Australian tables, the comparison 
revealed the following errors, all in S: 


x For Read 
0.05 .00298 12 .00297 30 
0.07 .00492 43 .00492 40 
1.97 .55507 4 .55507 3 


These three corrections were verified by the reviewer by independent compu- 
tation ; the error in the third case is about 1.3 final units. The compiler states that 
some values for small arguments were specially computed; it appears that this 
part of the work was not well done. On the other hand, it is surprising that a 6D 
table produced by subtabulating a 6D table should contain so few rounding errors 
in the portion tested ; perhaps the compiler knows more about the seventh decimal 
than he has claimed. The second differences appear to be modified when necessary, 
though this fact is not stated. The table does not pretend to be definitive, but will 
nevertheless be found very useful. 

A. F. 

1. G. N. Watson, A Treatise on the Theory of Bessel Functions, second edition, Cambridge 
a Press, England, 1944. 

2. S. Corrincton, Tables of Fresnel Integrals, Modified Fresnel Integrals, the Probability 


haa and Dawson's Integral, Radio Corporation of America, R.C.A. Victor Division. [MTAC, 
v. 7, UMT 166, 1953, p. 189.] 


1 = 
88[L ].—M. Ferentz & C. Harrison, “A tabulation of the function * f Jo(y)dy, 
0 
x = 0(.01)31; 4D,” Argonne National Laboratory, Lemont, Illinois, 8 ozalid 
sheets, 28 cm. Three copies deposited in the UMT FILE. 
This is a table of the function f(x) given in the title for x = 0(.01)31, to 4D. 


The table was computed from Taylor’s expansion. Values of f(x) are given for 


x = 0(.02)1 to 7D by G. N. Watson [1], p. 752. For x = 0(.02)16, 7D values of 
f(x) can be computed from tables given by Watson, l.c., p. 666 et seq., using 
the formula 


f(x) = Jo(x) + 5 5 (Ji(@)Ho() — Jo(x)H,(x)). 
P. HENRICI 
University of California 
Los Nee California 


G. N. Watson, A Treatise on the Theory of Bessel Functions, 2nd ed., Cambridge University 
Press, Cambridge, 1944. 
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s9[P].—R. E. D. Bishop & D. C. Jonnson, Vibration Analysis Tables, Cam- 
bridge University Press, New York, 1956, viii + 59 p., 28 cm. Price $2.00. 


In writing their book, Mechanics of Vibration, Cambridge University Press, 
the authors collected works from various previous engineering publications and 
managed, in part, the computation of formulae and numerical tables intended to 
help the analysis (and synthesis) of conservative mechanical systems. These aids, 
isolated from the main volume, have an independent value for a professional 
analyst; accordingly they have been extracted in this booklet form. 

The structural elements, of which the more complicated systems are assumed 
to be composed, and their motions are: lateral vibration of a taut, uniform string; 
torsional vibration of a uniform circular shaft ; longitudinal vibration of a uniform 
bar; flexural vibration of a uniform beam. For the three first cases, where the 
motion is governed by the same second order differential equation, the recep- 
tances, natural frequencies, and modes of vibration are given in terms of structural 
dimensions for three combinations of tiie simplest boundary conditions, namely, 
either the distortion or its longitudinal derivative vanishing at the ends. In giving 
the receptance formulae for the flexural vibration, the boundary conditions con- 
sidered are those of a beam with clamped, pinned, sliding, or free ends in all proper 
combinations. The related functions, composed of products and product sums of 
trigonometric and hyperbolic functions are given to 5S for x = 0(.05)11. For the 
boundary condition combinations, pinned-pinned, clamped-clamped, free-free, 
clamped-free, clamped-pinned, and free-pinned the characteristic functions and 
their derivatives up to the fifth mode are given to 5S for x/] = 0(0.02)1, and 
likewise the five first roots of the associated characteristic equations. 

These lucid tables form a valuable contribution to the similar and previous 
German and Russian collections. 


P. LAASONEN 
University of California 


Los Angeles, California 


90[U ].—Ernar ANDERSON (Director), TORBEN Krarup, & BJARNER SVEJGAARD, 
Geodetic Tables, International Ellipsoid, Geodaet. Inst. Skr. (Mémoires de 
l'Institut Géodésique de Danemark), Ser. 3, v. 24; Copenhagen, 1956, 8 p. 
introduction + 184 p. tables, 28.5 cm. 


Geodesists employ tables based on a standard assumed ellipsoid of revolution 
to extrapolate latitudes and longitudes from an astronomically determined point 
to new points. This procedure was originally worked out in the eighteenth century 
as a method of determining the figure of the earth by comparing a series of astro- 
nomical determinations with the results to be expected from some assumed shape 
of the earth. It was extended to the general mapping of the European countries 
as the original arcs of the meridian were expanded into the modern national 
triangulation nets. 

The present tables are for the most part carried to a precision of 12 decimal 
places or 12 significant figures. Since the question of the requirements for such 
precision have in the past been discussed in Mathematical Tables and Other Aids 
to Computation, it may not be out of place to discuss the problem here. 

The results of precise triangulation have a precision of the order of 073 in the 
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angles, or six-figure accuracy. Since this accuracy may be demanded between 
points only a few kilometers apart, it may, under some circumstances, imply 
accuracies of a few millimeters. Since 1 millimeter is 0700004 of latitude, the use 
of five decimals of seconds of arc may, under some circumstances, be required. 

At the same time, the triangulations of rather large areas, such as Europe or 
North America, have been brought to consistency by the simultaneous adjustment 
of the measured angles in very large systems of equations. Hence the latitudes and 
longitudes within the scheme must be treated formally as though they represented 
measurements with a precision of 10 or 11 significant figures. 

The problem comes to a head when it is necessary to transform the latitudes 
and longitudes into a coordinate system of the type which is used for military or 
cadastral surveys. There are only two ways to avoid the use of a large number of 
significant figures: first, to permit the precision of the computations to suffer, 
by rounding off. This loses the power of the coordinates to define the direction and 
distance between nearby points. Second, to make use of small systems, thereby 
risking trouble at the junctions. Neither is desirable; and hence the requirement 
for precise tables. 

It might be thought that the tables, although formally precise, need not be 
theoretically correct. An outstanding example of this scheme of thought is the 
tables for the Lambert Nord de Guerre, prepared in 1916 for France. Of the five 
constants at the head of this table, no three could be made to agree; and no two 
agreed with the table itself. This was not serious until, in 1943, it became neces- 
sary to extend the tables. One Allied agency extended them by differencing the 
tables and extrapolating ; another solved from the tables for the constants of the 
spheroid which would produce the given tables; while a third managed to guess 
the approximation which the originators of the table had made. Any one of these 
tables would have done the job, but the discrepancies between them meant that 
they could not be used together. Ultimately one was adopted. 

The present, very precise tables are based on the International Ellipsoid, 
adopted at Madrid in 1924, with the dimensions: 


a (equatorial radius) 6,378,388 meters 
f (flattening) 1/297. 


They replace logarithmic tables prepared at that time. They include: 
W = Vi — & sin? ¢, ¢ = 0°(1’)90°, 12D, 
107/N = 10’W/a, in meters; ¢ = 0°(1')90°, 12S, 
107/M = 10’w*/a(1 — e*) in meters“; ¢ = 0°(1’)90°, 12S, 


waC sin? 
tive | [ w + g.W 


@ — ¢*, defined by tan (¢*/2 + 2/4) = tan (¢/2 + 2/4) ( 


and expressed in seconds of arc; ¢ = 0°(1’)90°, 6D. 
8 — , defined by tan B = b/a tan ¢; b = aV1 — é 
and expressed in seconds of arc; ¢ = 0°(1')90°, 6D. 
p’/2MN, where p” = 648,000/x, ¢ = 0°(1')90°, 6S. 


| in milligals, ¢ = 0°(1’)90°, 8S, 


1 — esin +)" 
1+ esing 
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The quantities M and WN are the radii of curvature of the ellipsoid along the 
meridian and parallel to it. They are formed from the auxiliary quantity W. The 
terrestrial gravity 7, is here defined by the International gravity formula, slightly 
modified to make it rigorously consistent with the International ellipsoid. Note 
that the geodesists include centrifugal force in gravity, which they distinguish 
from gravitation, which is taken as not including centrifugal force. 

The Gaussian latitude, ¢*, is identical with the quantity called the isometric 
latitude in the tables of the U. S. Coast and Geodetic Survey and the U. S. Lake 
Survey. The term isometric latitude is used by the Europeans for gd-'¢*, the 
antigudermannian of the Gaussian latitude. The antigudermannian of ¢* is 
proportional to the north distance corresponding to the given latitude on a 
Mercator map. For a sphere, it reduces to the log tangent of the semi-colatitude. 
The Gaussian latitude is fundamental for the calculation of conformal projections 
(especially for coordinate computations). 

The reduced latitude, 8, is also called the parametric latitude. The meridional 
arc can be expressed as an elliptic integral of the second kind in terms of the 
parametric latitude. 

The tables have been computed from Fourier series by an IBM 602A, and 
printed directly from an IBM 416 tabulation. The figures are clear and readable. 

A cursory check of the table of ¢ — ¢* against the tables of latitude functions 
for the International Ellipsoid from the Lake Survey [1] indicated agreement to 
the 4th place of decimals, the limit of the Lake Survey tables. 


J. A. O'KEEFE 
Army Map Service 
Washington, D. C. 


1. War DEPARTMENT, Corps OF ENGINEERS, U. S. LAKE SuRvEy, Latitude Transformation, 
Hayford Spheroid, Geodetic Latitude to Isometric Latitude and Isometric Latitude to Geodetic Latitude, 
New York Office, Military Grid Unit, 1944. 


91[V].—SpPEER Propucts Company, The SPEER BALLISTICS CALCU- 
LATOR, Speer Products Company, Lewiston, Idaho, 23 X 10 cm. Price $1.00. 


This little instrument is a slide rule for the computation of drop and remaining 
velocity of small arms bullets in terms of muzzle velocity, range, and ballistic 
coefficient. The scales were prepared on the basis of Ingall’s tables, which were 
in turn based upon the Siacci approximation to flat trajectories and the GAvre 
drag function. As an example, the range scale (versus drop) goes from 50 to 1000 
yards, and may be read within about ten percent. Ballisti¢ coefficients are sup- 
plied for each of the manufacturer’s bullets. 

J. W. GREEN 


University of California 
Los Angeles, California 


92[V_].—HomeEr S. Pow.ey, Extension of Ingalls’ Table IIA. One typewritten 
sheet, 28 cm., deposited in the UMT Fixe. 


For Z = 20000(1000)45000 and 50000, the functions A V?/700?, log V/u and 
N are listed to 2D, 4D, and 2D respectively. For values of Z divisible by 5000 
A V?/700 is listed to 4D and N to 3D. 
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Ingall’s Tables are described by Bliss [1]. According to the author the present 
tables are intended to extend the tables for solution of problems relating to some 
low power weapons. 

ec B®. 


1. G. A. Biiss, Mathematics for Exterior Ballistics, John Wiley and Sons, New York, 1944, 
p. 33-35. 


93[W, X].—S. Vaypa, The Theory of Games and Linear Programming, John Wiley 
& Sons, Inc., New York, 1956, 106 p., 17 cm. Price $1.75. 


This handy pocket-sized book contains an exposition of some of the most useful 
aspects of the theory of two-person zero-sum games and of linear programming. 

Any work of this size must necessarily be incomplete. Thus the present volume 
contains no reports of experience with extensive computations involved with 
games and linear programming problems. There is no indication of the handling 
of non-linear situations, and no real full use of the continuous space in which the 
problems are invented. A final chapter devoted to Beale’s “‘Method of Leading 
Variables,” and a few introductory and incomplete historical remarks miglit 
profitably have been omitted in favor of other topics (the genesis of problems, for 
example), but on the whole there can be no real objections to the choice of 
material expounded. 

A major effort is devoted to a simplex method, which is certainly the most 
widely used method resolving games and linear programming problems. The 
description of the theory covered is elementary and clear. However there are no 
exhibited statements, such as theorems, to summarize the arguments that come 
to mind. 

Several examples are given and solved. The number of drawings included to 
clarify the presentation is impressive, particularly in the early chapters which are 
devoted to graphical representations of the theory. This all contributes to the 
exposition. 

In all the book is a valuable contribution to the literature in this popular field, 
and it contains a valuable selection of material presented clearly. 

ce ae ee 


94[P, W, X, Z ].—Proceedings of the Second Annual Computer Applications Sym- 
posium, held October 24-25, 1955, sponsored by the Armour Research Founda- 
tion of the Illinois Institute of Technology, Chicago, Illinois, 1956, 108 p., 
23 cm. Price $3.00. 


This is a collection of reports of talks delivered at the symposium, together 
with the discussion which followed each paper. For all practical purposes, the 
word ‘‘digital” could be in the title of the symposium. 

The program was devoted to seven papers having to do with ‘Computers for 
Business and Management,” and seven concerned with ‘‘Computers for Engineer- 
ing and Research.”” Twelve papers and one abstract appear in the Proceedings, 
as follows: 

Computers for Business and Management: ‘The use of digital computers in 
industry,” by R. F. Clippinger, ‘A dollar and cents approach to electronics,”’ by 
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John L. Marley, “An application of computers to general bookkeeping,” by 
W. F. Otterstrom, “User experiences and applications of the ERA 1103,” by 
George E. Clark, ‘“‘Automobile selective underwriting and automatic rating on 
the IBM 650,” by C. A. Marquardt, ‘‘Cutting costs with linear programming,” 
by Jacob E. Bearman, ‘‘Probability forecasts in management decisions (abstract),” 
by Stanley Reiter. 

Computers for Engineering and Research: ‘‘Use of the IBM 650 in scientific 
computations,” by A. W. Wymore, “Engineering applications of large scale com- 
puters,” by C. B. Ludwig, “High speed computation of engine performance,” by 
J. T. Horner, ‘Pyrolysis reactor design computations,’ by H. C. Schutt and 
R. H. Snow, “Aircraft flight test data processing,’’ by T. M. Bellan, ‘“‘Program- 
ming a Monte Carlo problem,”’ by J. F. Hall and J. M. Cook. 

The Proceedings contain some information which is very useful for those 
interested in the two broad fields covered. 

T. H. SouTHarp 


University of California 
Los Angeles, California 


95(G, H, X].—Ruicuarp B. Smitu, Table of Inverses of Two IIl-Conditioned 
Matrices, Westinghouse Electric Corporation, Bettis Atomic Power Division, 
Pittsburgh, Pennsylvania, 1957, ii + 68 p., 29cm. Deposited in the UMT Fine. 


The listings are the inverses (m = 2,3, ---, 15) for A, = (aij), a1; = 1 for 
j = 1,2, ---,n, a5 = 1/(¢ +7 — 1) fori = 2,3, ---,m andj = 1, 2, ---, m; and 
B, = (bij), bs5 = 1/(6 +72+ 7 — 1) for p = Oandi,7 = 1, 2, ---, m. The matrix 
A, is described by M. Lotkin [1] and B, is described by Savage and Lukacs [2] 
and Collar [3, 4]. 

The calculations were carried out on an IBM 650 using a complete double 
precision rational number interpretive system. The listings are believed to be 
correct. 


RICHARD B. SMITH 
Westinghouse Electric Corporation 
Bettis Atomic Power Division 
Pittsburgh, Pennsylvania 

1. Mark Lorn, “‘A set of test matrices,” MTAC, v. 9, 1955, p. 153-161. 

2. I. R. Savace & E. Luxacs, “Tables of inverses of finite segments of the Hilbert matrix,” 
NBS Applied Mathematics Series, No. 39, Contributions to the Solution of Systems of Linear Equa- 
tions and the Determination of Eigenvalues, 1954, p. 105-108. 

3. A. R. Cotzar, “On the reciprocation of certain matrices,’’ Royal Soc. Edinburgh, Proc., 
v. oy we p. 195-206. 

A.R. Co.xar, “On the reciprocal of a segment of a generalized Hilbert matrix,’”’ Cambridge 
Phil “Soc., Proc., v. 47, 1951, p. 11-17. 


96[H, X].—Zycmunt Doweirp, Krakowiany i ich zastosowanie w mechanice 
budowli. (Cracovians and their application in structural mechanics.) Pahstwowe 
Wydawnictwo Naukowe, Warszawa, 1956, 168 p. zt. 18. 


Cracovians are rectangular arrays of numbers which are added like matrices, 
but which are multiplied column-by-column. They were developed by T. 
Banachiewicz, apparently because column-by-column multiplication is easier in 
desk computing than row-by-column multiplication. The present work expounds 
the definitions, notations, and properties of cracovian theory, and their use in 
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problems of linear algebra. The second chapter is devoted to the solution of sys- 
tems of linear equations by various methods of triangular decomposition. In the 
third chapter are discussed simple iterative methods for solving linear systems 
and for computing eigenvalues. The exposition is elementary. There are problems 
from structural mechanics, and many numerical examples of orders up to 5 or 6. 
The examples are oriented towards desk computation. 

The fact that cracovian multiplication is non-associative causes various 
strained notations, and appears to the reviewer as an overwhelming impediment 
to fundamental progress. Nevertheless, cracovians have a minority of enthusiastic 
supporters, largely but not exclusively in Poland. 


G. E. ForsyTHE 
University of California 


Los Angeles, California 


(After September, 1957, will be at Stanford University, Stanford, California.) 
This review was prepared by G. E. Forsythe for Mathematical Reviews. 


97[S, X ].—Srr HAROLD JEFFREYS & BERTHA SwIRLEs (LaDy JEFFREYS), Methods 
of Mathematical Physics, Cambridge University Press, Great Britain, 1956, 
iv + 714 p., 25 cm. Price $15.00. 


This is the third improved edition of this impressive textbook of mathematics 
as it should be applied to physics. The authors are firm in their feeling that the 
physicist needs a rigorous proof of all theorems used, and they present such proofs 
under conditions which are suitable for the physics applications they have in mind. 

At the beginning of the chapter on numerical methods the authors quote 
Lord Kelvin, ‘‘I have no satisfaction in formulas unless I feel their numerical 
magnitude.” For purposes of this book the authors have no interest in mathe- 
matical theorems unless they are applicable to physics (their stated standard is 
applicability in at least two branches), and they have no satisfaction unless the 
theorem is proved under adequate conditions for their application and an appli- 
cation illustrated. 

Numerous problems are included. 

Much of the work concerns numerical analysis; this varies from the classical 
studies of finite differences through various uses of analysis and studies of special 
functions. Thus the book is a valuable but not a complete textbook on many 
aspects of numerical analysis. 

The authors have continued to improve the book through these three editions 
evidently seeking advice wherever it is available. Thus they have put together 
a sound book containing much material not easily available elsewhere. 

The chapter headings follow: 

Chapter 1. The Real Variable, 2. Scalars and Vectors, 3. Tensors, 4. Matrices, 
5. Multiple Integrals, 6. Potential Theory, 7. Operational Methods, 8. Physical 
Applications of the Operational Method, 9. Numerical Methods, 10. Calculus of 
Variations, 11. Functions of a Complex Variable, 12. Contour Integration and 
Bromwich’s Integral, 13. Conformal Representation, 14. Fourier’s Theorem, 15. 
The Factorial and Related Functions, 16. Solution of Linear Differential Equa- 
tions of the Second Order, 17. Asymptotic Expansions, 18. The Equations of 
Potential, Waves, and Heat Conduction, 19. Waves in One Dimension and Waves 
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with Spherical Symmetry, 20. Conduction of Heat in One and Three Dimensions, 
21. Bessel Functions, 22. Applications of Bessel Functions, 23. The Confluent 
Hypergeometric Function, 24. Legendre Functions and Associated Functions, 
25. Elliptic Functions. 
There is a reasonably detailed index. 
C. B. T. 


98[S, X].—F. B. HILDEBRAND, Advanced Calculus for Engineers, Prentice-Hall, 
Inc., New Jersey, 1956, xiii + 594 p., 22 cm. Price $7.75. 


The reviewer belongs to a school of thought which holds that the simplest 
numerical methods (Euler’s method of forward differences in ordinary differential 
equations, for example) may be the soundest introduction to many mathematical 
subjects. He would have liked to find more numerical material in this text by 
this author of one of our better numerical analysis texts. 

The volume does contain the following tables: A table of 41 Laplace transforms, 


a table of RF x"Jm(x) and vf x"Im(x) form = —4}(1)9/2, P(x), x = 1(.01)1.99, 


4D, and the first five zeros of J,(x) p = 0(1)5. 

A chapter on numerical methods for solving ordinary differential equations 
includes Taylor’s series, the Adams method, the Runge-Kutta method, and the 
Picard method. A chapter devoted to series solutions of differential equations 
introduces Bessel functions, Legendre functions, and the hypergeometric function 
(and the gamma is introduced in connection with Laplace transforms). 

Various iterative methods for computing eigenvalues or solutions of equations 
are discussed, and the author does pay attention to the requirement for numerical 
answers which engineers frequently face. However, when treating partial differ- 
ential equations the author does not take up difference equation methods of 
approximating solutions, which were not popular at the time of writing. 

Several topics usually included in advanced calculus texts are omitted; 
multiple integrals and surface integrals are treated only as they arise in vector 
analysis. 

Problems are carefully chosen. 

The book more or less parallels material in the much more ambitious and 
difficult work on Methods of Mathematical Physics, by Jeffreys and Jeffreys [1]. 
Chapter headings are listed below. 


Solutions of Linear Ordinary Differential Equations , 
The Laplace Transformation 
Numerical Methods for Solving Ordinary Differential Equations 
Series Solutions of Differential Equations 
Boundary-Value Problems and Orthogonal Functions 
Vector Analysis 
Partial Differential Equations 
Solutions of Partial Differential Equations of Mathematical Physics 
Functions of a Complex Variable. 
oy Hy 


i. Maaee> jerome & BERTHA SWIRLES JEFFREYS, Methods of Mathematical Physics, Cam- 
bridge, England, at the University Press, 1946. 
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99[X, P].—S. H. CRANDALL, Engineering Analysis, A Survey of Numerical Pro- 
cedures, McGraw-Hill Book Co., Inc., New York, 1956, x + 417 p., 24 cm. 
Price $9.50. 


This textbook has grown out of the author’s continuing efforts to acquaint 
engineering students with numerical methods of attacking problems in engineering 
analysis. The author describes engineering analysis as the performance of two 
steps: ‘‘1. Construction of a mathematical model for a physical situation. 2. Re- 
duction of the mathematical problem to a numerical procedure.” 

The present book is largely devoted to the second of these steps, but the ar- 
rangement of material conforms more nearly with the first step than with the 
arrangement of mathematics courses in most schools. Actually, to a mathe- 
matician, the book is more nearly a most valuable catalogue of methods (with 
careful references to existing literature) than a textbook, for proofs are presented 
only in summary form. 

Three classes of problems are considered: Equilibrium problems, Eigenvalue 
problems, and Propagation problems. Each of these classes is considered first as 
a lumped-parameter problem (or a problem with a finite number of variables) 
and then (three chapters later) as a continuous problem. 

The reviewer has long felt that use of at least the most elementary numerical 
methods is the soundest introduction to many courses in mathematics. Euler's 
method of forward differences for numerical solution of an ordinary differential 
equation, for example, offers the student a feeling for the properties described in 
the Picard existence and uniqueness theorem ; calculation of difference quotients 
(particularly when they are stated in terms of displacement and time) seems to 
be a sound introduction to the derivative, and so on. Thus one important purpose 
which this book can serve is presenting this feeling concerning the nature of 
solutions of the problems treated. Just as a student who knows Euler’s method 
can do something with any ordinary differential equation initial value problem 
he is likely to meet, so the student of the present text may make efficient progress 
on any problem he is likely to encounter from the fields studied. The study (or 
at least the perusal) of a book of this type would seem to provide an excellent 
introduction to many aspects of abstract analysis through the concrete examples 
furnished. Here the motivation for the non-numerical studies would be the promise 
that some of the problems can be solved more generally and with less effort—a 
situation to which many students seem highly attracted. 

In any event, since recent engineering problems have more and more demanded 
numerical solution, the book (or its equivalent, which seems not to exist) seems 
necessary for any complete training in several fields of engineering, and the re- 
viewer suggests that teachers of more abstract courses might well extract much 
material from the book for introductory use. 

The author has chosen methods which he feels are most likely to be useful 
in engineering analysis, and he’ has expounded them carefully and in a scholarly 
way consistent with the level of maturity at which he aims. In approximate 
methods he mentions stability considerations carefully, but he does not swamp 
the student with all technicalities which can crop up in stability studies. In con- 
nection with relaxation he mentions overrelaxation, but stops short of the detailed 
studies available on the subject. 
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The material to be covered requires a study of methods of solution of systems 
of linear algebraic equations, and ordinary and partial differential equations with 
various kinds of initial and boundary conditions. Numerical methods (of scope 
indicated below) are presented in detail along with an outline of the basic mathe- 
matics involved. 

The material in the book is illustrated with many worked examples and with 
a spectacular number of figures and much tabulated material. Problems for the 
student in each section vary from routine application to problems demanding 
some knowledge of the theory—determination of degree of convergence of discrete 
methods, for example. 

There is an impressive number of references to basic material used in the book 
and to additional material which is available for more extensive study. 

The book is the second which we have seen recently stressing the importance 
of numerical calculation in engineering analysis; Purday’s book [1] covers much 
of the same material but with considerably less detail. 

A reasonable idea of the contents is gained from the section headings, which 
follow. 


1. Equilibrium problems in systems with a finite number of degrees of freedom: 
Particular examples—Formulation of the general problem—Mathematical prop- 
erties—Extremum problems—Elimination method for linear systems—Iteration 
—Relaxation—Iteration combined with elimination—Procedures applicable 
directly to physical systems. 

2. Eigenvalue problems for systems with a finite number of degrees of freedom: 
Particular examples—Matrix notation—Formulation of the general problem— 
Mathematical properties—An extremum principle for eigenvalues—Direct 
methods of solution—Iteration—Intermediate eigenvalues—n-step iteration— 
Relaxation methods—Upper and lower bounds for eigenvalues—Diagonalization 
of matrices by successive rotations. 

3. Propagation problems in systems with a finite number of degrees of freedom: 
Particular examples—Formulation of the general problem—Mathematical prop- 
erties—Iteration—Series methods—Trial solutions with undetermined parameters 
—Finite-increment techniques—Introduction to step-by-step integration pro- 
cedures—Recurrence formulas with higher-order truncation error—Step-by-step 
integration methods for systems with several degrees of freedom. 

4. Equilibrium problems in continuous systems: Particulax examples—Formula- 
tion of the general problem—Mathematical properties—Extremum problems— 
Trial solutions with undetermined parameters—Finite-difference methods. 

5. Eigenvalue problems in continuous systems: Particular examples—Formula- 
tion of the general problem—Mathematical properties—Extremum principles for 
eigenvalues—Iteration—Trial solutions with undetermined parameters—Finite- 
difference methods. 

6. Propagation problems in continuous systems: Particular examples—Formula- 
tion of the general problem—Mathematical properties—Trial solutions with un- 
determined parameters—Finite-difference methods for parabolic systems—Finite- 
difference methods for hyperbolic systems. 





of | 
eng 
to s 
tior 
to t 


tio 


ch: 


nic 





stems 
3 with 
scope 
1athe- 


| with 
or the 
nding 
screte 


- book 


‘tance 
much 


which 


edom: 
prop- 
ration 
icable 


edom: 
lem— 
Direct 
tion— 
zation 


Edom: 

prop- 
neters 
 pro- 
y-step 


‘mula- 
ems— 


‘mula- 
les for 
*inite- 


‘mula- 
th un- 
*inite- 





REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 221 


One note of criticism concerns the first sentence of the preface: ‘“‘The advent 
of high-speed automatic computing machines is making possible the solution of 
engineering problems of great complexity.’”’ However, references in the book 
to such machines or experience with machines are sparse. This valuable publica- 
tion certainly did not need this implicit and unfulfilled promise of an introduction 
to the bright new electronic world. 

Ce: &. 


F. P. Purpay, Linear Equations in Applied Mechanics, Interscience Publishers, Inc., 
New “York, 1954. [Rev. 66, MTAC, v. 9, p. 131-4 


100[X, Z ].—GeorceE R. Strsitz & JuLes A. LARRIVEE, Mathematics and Com- 
puters, McGraw-Hill Book Co., Inc., New York, 1957, vi + 228 p., 23 cm. 
Price $5.00. 


The eleven chapters of this book cover nearly all phases of digital computer 
design and use. The first three chapters cover the basic elements of computers 
and mathematics, such as the basic differences between the analog and digital 
computers and the mathematical definition of a function. These early chapters 
also include some of the history of mathematics and computing, and outline the 
types of problems arising in applied mathematics which can be and are solved 
by computing techniques. The fourth chapter is a discussion of the history of 
computers which is both interesting and informative but, unfortunately, does 
not go beyond the ENIAC computer to a description and discussion of the more 
recent history of stored program computers. The fifth chapter contains discussions 
of some numerical techniques including those for finding roots of polynomials, 
solving linear equations, and solving ordinary and partial differential equations. 
This chapter gives an excellent insight to the beginning student on some of the 
techniques, but many of the techniques presented are not practical for the modern 
electronic computer. Bernoulli’s method for finding roots of a polynomial equa- 
tion, for example, is of little importance in the modern computing world. The 
next three chapters discuss digital computer components, number systems, com- 
puter memories, and analog-digital converters. These chapters also outline tech- 
niques for performing arithmetic processes in the digital computer, and contain 
an elementary description of how the computer carries out its operations on the 
basis of a stored program. The analog computer is likewise treated in these 
chapters, and the techniques for the mechanical differential analyzer in performing 
integration are explained. Chapter 9 is a discussion of Monte Carlo techniques 
and includes an explanation of the solution of linear equations by these tech- 
niques. The discussion on sampling techniques includes brief discussions on the 
solution of certain partial differential equations and techniques for generating 
“random’’ sequences. Chapter 10 is a discussion on the various errors which can 
occur in both analog and digital computers. Chapter 11 summarizes various 
special-purpose applications of computers such as those seen in the automatic 
factory, language translation, and, in a lighter vein, computers which play games. 

The central difficulty of this reviewer in evaluating this book is in deciding 
to whom the information is directed. The book seems to be too elementary to 
serve as a reference book for any computer professional, either user or designer. 
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On the other hand, the book does not go into sufficient detail in any area to be a 
successful textbook for most university courses. It could, however, serve as a 
textbook on an introduction to computers given to freshman or sophomore stu- 
dents or more advanced non-science students. It could also serve as a reference 
for the beginner in the field or for a non-professional. 

The most serious shortcoming of the book is that the material is not sufficiently 
modern. Many of the numerical techniques presented are old-fashioned and funda- 
mentals are often described in terms of old-fashioned equipment rather than 
modern equipment. (An analog-digital converter could be easily described in 
terms of halving and comparing the voltage rather than in terms of the mechanical 
device which transmits a shaft position to a mechanical digit position.) The 
authors leave the false impression that tables, especially those involving the 
elementary functions, are frequently stored in computer memory. Often modern 
words are not used; for example, the words “subroutine,” “‘programmer,” and 
the term “parity check bit,” are not used in describing and discussing these items. 

The style of the book is light and enjoyable and makes interesting reading. 
The bibliography is good and extensive. 

WALTER F. BAUER 


The Ramo-Wooldridge Corp. 
Los Angeles, California 


101[X, Z].—W. J. Eckert & REBEccA JonEs, Faster, Faster, McGraw-Hill Book 

Co. Inc., New York, 1955, vii + 160 p., 23 cm. Price $3.75. 

Quoting from the preface, “This monograph is an attempt to explain in non- 
technical language how a calculator operates, the nature of the problems it solves, 
and how the probiems are presented to the calculator.”’ Actually, it consists almost 
entirely of a description of ““NORC,” the Naval Ordnance Research Calculator 
designed and built by the International Business Machines Corporation. Probably 
but few experts would agree to the claim that NORC “. . . is also easiest to 
understand” (preface) yet the explanatory attempt must be judged very suc- 
cessful; it is completely devoid of engineering details (such as component types, 
circuit diagrams, etc.), relying instead upon simple block diagrams and schematics 
accompanied by descriptions of the essential properties involved. The language 
is indeed nontechnical, even such common descriptive contractions as ‘‘and/or 
gate” being avoided, although various colloquialisms native to business ma- 
chinery, such as “‘echo-pulse,”’ ‘‘print cycle,” are explained and used. Here and 
there the viewpoint tends to be a bit insular; for instance the basic electronic 
building block turns out to be a binary pulse shaper producing 1 microsecond time 
delay ; this is called a Dynamic Pulse circuit, written always with capitals, like 
a Thing. Again, no allusion is made to any except IBM equipment, nor to tech- 
nical contributions by any outsiders except F. C. Williams, the Red Queen, 
Leibniz and Newton. 

Chapter I is introductory, and begins with a general outline of the need for, 
and requirements of automatic en masse arithmetic, and sketches the typical 
mental, organizational, and instrumental hurdles that must be overcome. Lucid 
examples of basic arithmetical and scaling operations are given, as well as various 
notions about electronic components, timing relationships and the réle of main 
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functional units such as the arithmetical, memory, input-output, etc. Halfway 
through the chapter all such generality is jettisoned without rites and supplanted 
by a clear outline of NORC’S code-language, organization, structure, and 
philosophy ; these and other features of NORC are more fully developed in all 
chapters but the last of the ten that follow, and in five summarizing appendices. 
The level of exposition is, on the whole, excellent; one may take exception—as 
does this reviewer—to the practice of attaching grammatical commas and periods 
to arithmetical examples, and to the use of the word “accuracy” instead of the 
word “reliability’”’ in contexts where freedom from accidental errors is meant, 
but aside from these trivialities of taste, the exposition leaves little to be desired. 

NORC is revealed to be quite a fast machine (70 microseconds to multiply; 
50 to add or subtract ; 250 to divide; 8 for memory access), with ample memory 
capacity (3600 words, each 66 bits) composed of 4 banks of Williams CRT storage, 
supplemented by 8 magnetic tape units of impressive capacity (400,000 words 
each) and speed (4000 words/second), plus printing (19 words/sec) and card 
punch/read facilities (400 words/min). Each order is an instruction plus three 
addresses, the latter automatically modifiable en route. The instruction list is 
almost lavishly flexible; all standard arithmetic operations are available with 
options of floating point (+30), automatic sizing shifts, extractions for order 
modification, transplantation of addressees, etc. There are no less than 21 con- 
trol transfer instructions, plus 5 for print control and 9 for tape control; the 
grand total is 98. Data entries are normally preserved to 13 decimal places, and 
there appears to be an efficient provision for double-precision arithmetic. Beyond 
doubt, NORC should be classed as an outstandingly flexible and effective com- 
puting system. 

Electronically, NORC belongs in the category of “‘pulse recirculating’’ ma- 
chines typified by SEAC of the National Bureau of Standards. These machines 
avoid the use of static memory cells of Eccles-Jordan “flip-flop” type in their 
arithmetical and logical units by routing binary pulses through various circulating 
paths consisting of time delay and pulse re-shaping components. Since identifica- 
tion of pulse position on this moving coordinate system is vital, a central clock 
is used to quantize the time coordinate and great stress is laid upon matters of 
timing, pulse shaping, and synchronism. 

The internal language of NORC consists of words 16} decimal digits long, each 
digit being represented by a tetrad of binary digits. It is actually true that arith- 
metic within the machine is carried out in this mixed base, elements of an essen- 
tially binary nature (aggregates of dynamic pulse units plus gangswitches) being 
combined so as to preserve the local identity (i.e., within tetrads of wires, ele- 
ments, etc.) of decimal characters. Words are taken from or put into the memory 
as 66-bit parallel transfers, whereas within the arithmetic unit the basic process 
involves 4-abreast shifting of successive tetrads to effect, for instance, addition 
that is decimally serial. Multiplication is done by a similar shift of the multiplicand 
stepwise through a parallel circuit called a ‘product generator’ equivalent to a 
multiplication table, so that as each multiplicand digit steps up to the altar, all 
ten multiples of it become simultaneously available. From this array all multiples 
specified by the multiplier digits can thus be accumulated at each multiplicand 
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step, and the entire product assembled during one ‘‘pass.’’ Consequently multi- 
plication is almost as brief as addition, though division does not fare comparably 
well, requiring four or five times as long. 

A little bookkeeping may suggest that NORC is somewhat extravagant with 
regard to information capacity, one binary tetrad capable of representing 16 
alternatives being used to represent only 10, etc. This glaring redundancy is the 
decimal man’s burden and seems unavoidable; in NORC some slight byproduct 
utility is recovered by using ‘12’’ and “13” for magnetic tape word-end and 
block-end signals. Much emphasis is placed upon automatic checking throughout 
the discussion and NORC uses two systems, (a) modulo 4 summing of the bits 
in each word, affixing half a tetrad of redundancy to tag this, and (b) “‘casting- 
out-nines”’ checking of the arithmetic operations. Together, (a) and (b) seem to 
make good sense, whereas (a) alone would be quite weak since repeated doublings 
play a key réle in the product generator. In the Williams memory a parity check 
is also provided for the sum of the bits in each of the 66 parallel positions. Alto- 
gether, these features inspire confidence, yet in a machine having some 3 of its 
capacity redundant, one wonders whether it might not have been feasible to make 
the checking density far more severe. 

The reader concerned with machine design will appreciate that NORC has 
been moulded to fit rigorously within a framework of precepts: decimal, speedy, 
ample capacity, big vocabulary, pulsed circuits, checking. Probably not everybody 
would choose exactly this prescription, but few would deny that it is interesting 
and bold, and seems to have been carried out with systematic skill and generous 
material resources, and with the sort of opportunistic ingenuity that is the hall- 
mark of elegant design. The field of computer design is advancing so rapidly both 
in technology and in concepts that no design group can claim ultimate insight, 
and all may well benefit from a study of the various species evolved. NORC 
deserves careful study, and it is to be hoped that a careful comparative examina- 
tion of its operational performance will eventually be published. 

The final chapter in Faster, Faster is entitled, ‘‘What is there to Calculate?”, 
and consists of a very elementary but clear discussion of calculations applied to 
linear systems, ballistic trajectores, planetary motion, etc. Perhaps a slight im- 
provement could be made by distinguishing more carefully the idea of solving a 
mathematical problem in general from the idea of calculating a numerical solu- 
tion; on page 131 following a discussion of ballistics, the statement that “A 
slightly more complicated problem is the three-body problem”’ may illustrate this 
point. Aside from this, the discussion of numerical procedures is well suited to 
convey the flavor of this field to the nonspecialist. 

J. H. B. 


102(F, Z].—D. H. Lenmer, “Sorting cards with respect to a modulus,” J. Assn. 
for Comp. Machinery, v. 4, 1957, p. 41-46. 


The author gives a method of using an ordinary punched card sorter to 
separate a deck of punched cards into no more than m sub-decks in each of which 
all numbers are congruent mod m. 

Cenk 
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103[Z ].—L. Pease, An Improved Control Board for Card-Programming the IBM 
602A Calculating Punch, Atomic Energy of Canada, Limited, Chalk River, 
Ontario, A.E.C.L. Report No. 317, 1955, reprinted 1956, vi + 24 + 8 p. of 
diagrams and figures, 27 cm. Price $1.00. 


The paper gives detailed wiring diagrams and instructions for use of a general 
purpose 602A board. Although the 602A is an electromechanical and, by electronic 
standards, a slow machine, it has a high degree of flexibility which is utilized in 
the board described. The board will perform operations on 8-digit numbers of 
the form: (A + a)*(B + b) + c = d where * indicates +, —, X or /; the small 
letters indicate any of seven internal storage registers, and the capital letters 
fields on the input cards. 


OweEN R. Mock 
North American Aviation, Inc. 
Inglewood, California 


104[Z ].—W. G. BROMBACHER, JULIAN F. SmitH, & Lyman M. VAN DER PYL, 
Guide to Instrumentation Literature, NBS Circular 567, U. S. Gov. Printing 
Office, Washington, D. C., 1955, iv + 156 p., 26 cm. Price $1.00. 


A bibliography listing pertinent works in instrumentation. It includes sections 
on automatic control, analogue and digital computers, and many other applica- 
tions more or less related to computation. Books and reference works which 
appeared no more than about twenty years prior to the compilation of the bibli- 
ography were listed, and all periodical articles found were listed. The closing date 
is not stated, but it was presumably some time during 1954. Sources of material 
are generally described and listed as specifically as is reasonable. 

i OA 


105[W, Z].—NaTIONAL PuysicaL LaBoratory, Wage Accounting by Electronic 
Computer, Report No. 1 of the Inter-Departmental Study Group on the 
Application of Computer Techniques to Clerical Work. Her Majesty’s Sta- 
tionery Office, London, 1956, 25 cm. Price 2s. 6d. net. 


This 57-page booklet reports on one of the earliest applications of computers 
to commercial work in England; the application is a government payroll calcula- 
tion, and the machine used is the DEUCE. The report is quite complete, including 
block diagrams, time and cost data, and even considerable discussion of computers 
per se, reliability, and a section on data sorting. Following are some of the report’s 
conclusions : 


(1) A “scientific computer’’ with magnetic tape input-output forms an ade- 
quate machine for payroll work. 

(2) The reliability of this operation is satisfactory if attention is given to 
program checks, proper machine maintenance, and safety margins in the time 
schedule ; automatic ‘‘built-in’’ checking is not essential. 

(3) The computer permits reduction in clerical staff, but further economies 
could be obtained if input data were originally recorded in appropriate form. 

(4) The report contains cautious statements indicating that payroll work may 
justify or at least help justify a computer installation in many organizations. 
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The DEUCE (successor to ACE) has a 250,000 bit magnetic drum, 32 bit 

word size, mercury delay line high speed store, and two milliseconds multiply time. 
D. D. Watt 


IBM Corporation 
Los Angeles, California 


TABLE ERRATA 
The following errata are mentioned in this issue: 


CARL-ERIK FROBERG, Hexadecimal Conversion Tables, Review 82, p. 208. 

H. NAGLer, Table of Square Roots of Integers, Review 77, p. 205. 

T. Pearcey, Table of the Fresnel Integral to Six Decimal Places, Review 87, 
p. 210-211. 

J. RyBNER, Nomogrammer over komplekse hyperbolski funktioner, Review 80, 
p. 207. 

G. N. Watson, A Treatise on the Theory of Bessel Functions (1. M. Longman 
paper, p. 179]. 


256.—GEORGE WELLINGTON SPENCELEY, RHEBA MuRRAY SPENCELEY, & EUGENE 
RHODES EPPERSON, Smithsonian Logarithmic Tables to Base e and Base 10, 
The Smithsonian Institution, Washington, D. C., 1952. [Review 992, MTAC, 
v. 6, 1952, p. 150-151. ] 


On p. 241 for log 1902 
read log 1902 


3,27921 05129 01395 12706 
3,27921 05126 01395 12706. 
; J. RAFALOwICcz 
B. JAKUBOWSKI 


Il 


Dept. of Physics 
Technical Institute of Wroclawska 
Wroclawska, Poland 


NOTES 
Handbook of Mathematical Tables 
National Bureau of Standards 


The National Science Foundation has commissioned the National Bureau of 
Standards Applied Mathematics Division to prepare a Handbook of Mathematical 
Tables containing formulas and graphs. This preject is an outgrowth of a con- 
ference on Mathematical Tables held at Massachusetts Institute of Technology 
on September 15 and 16, 1954. One of the principal recommendations made at 
this conference was that “‘an outstanding need is for a ‘Computer’s Handbook,’ 
with usually encountered functions, together with a discussion of their analytic 
properties and a set of formulas and tables for interpolation and other techniques 
useful to the occasional computer.” 
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Subsequently Dr. P. A. Smith, Chairman of the Mathematics Department of 
the National Research Council, appointed a committee composed of A. Erdélyi, 
M. C. Gray, N. C. Metropolis, P. M. Morse (Chairman), R. D. Richtmeyer, 
J. B. Rosser, H. C. Thacher, Jr., John Todd, C. B. Tompkins, and J. W. Tukey 
to advise the NBS staff and help establish a philosophy for the Handbook. After 
several meetings the following points were established : 


1. The proposed preliminary table of contents is that there shall be a general 
introductory section explaining the use of the tables and the following chapters: 
1. Powers, Roots and Related Functions; 2. Binomial Coefficients, Bernoulli Numbers; 
3. Fundamental Constants; 4. Circular and Hyperbolic Functions, arithms; 5. Sine, 
Cosine, Exponential and Logarithmic Integrals; 6. Gamma and Related Functions; 7. Error 
Function, Fresnel Integral ; 8. Legendre Functions; 9. Bessel Functions including: a) Integral 
Order, b) Spherical, Modified Spherical, Fractional Order, c) Complex ment, d) Integrals, 
e) Struve Functions; 10. Elliptic Functions and Integrals; 11. Mathieu Functions, Spheroidal 
Wave Functions; 12. Parabolic Cylinder Functions; 13. Hypergeometric Functions; 14. 
Confluent Hypergeometric Functions; 15. Miscellaneous Functions; 16. Orthogonal Poly- 
nomials; 17. Statistical Tables; 18. Interpolation Coefficients—Quadrature; 19. Radix 
Conversion Tables; 20. Combinatorial Tables. 


2. The tables are not to be given to a uniform number of decimal places or 
significant figures. The elementary functions are to be given to a high order of 
accuracy because they have basic importance and are essential for the evaluation 
of higher mathematical functions with the aid of auxiliary functions. 

3. The tables will be linearly interpolable to 5D or 5S, as far as possible, 
although more figures may be given. 

4. Polynomial and rational approximations for the various functions are to 
be given as interpolation and computing aids. 

5. Graphs are to be given to demonstrate the behavior of the function or as 
means of tabulation in the case of some of the higher transcendental functions. 

6. Auxiliary functions and arguments will be used extensively to simplify 
interpolation and permit tabulation over the entire range of the argument. 

7. Radix and “Key values’’ tables to high order of accuracy will be given. 

8. In each of the chapters there will be given the most useful mathematical 
relations pertinent to the functions, indefinite and definite integrals, infinite series, 
inequalities, and references supplementing the tables. 


Part of the computations have already been completed and it is hoped that 
the volume will be completed in about fifteen months. Dr. M. Abramowitz, of 
the Applied Mathematics Division of the Bureau of Standards, is in charge of 
the preparation of the Tables. Queries and suggestions should be addressed to him. 

Puitip M. Morse 


Massachusetts Institute of Technology 
Cambridge, Massachusetts 


International Congress of Mathematicians 
August 14-21, 1958 


The International Congress of Mathematicians will meet in Edinburgh, 
Scotland from August 14th to August 21st, 1958. The Executive Committee is 
inviting a number of mathematicians to deliver one-hour and half-hour addresses. 
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There will also be daily sessions devoted to fifteen-minute communications. Th 
will be eight sections, namely : 


. Logic and Foundations 

. Algebra and Theory of Numbers 

. Analysis 

. Topology 

. Geometry 

. Probability and Statistics 

. Applied Mathematics, Mathematical Physics, and Numerical Analysis 
8. History and Education 


Those who wish to receive further information about the Congress are 
quested to communicate their names and full addresses to the Secretary, Franl 
Smithies at the Mathematical Institute, 16 Chambers Street, Edinburgh 
Scotland. 

(Extracted from preliminary announcement) 


N.P.L. Mathematical Table Series 


The National Physical Laboratory has announced an N.P.L. Mathematie 
Table Series to contain tables arising from computational problems received i 
the Mathematics Division of the N.P.L. 

The first volume, The Use and Construction of Mathematical Tables, by L. Fox 
will be reviewed shortly in MTAC. 

Forthcoming volumes which have been announced are: 


Vol. 2, L. Fox, Tables of Everett Interpolation Coefficients, and 
Vol. 3, G. F. Miller, Tables of Generalized Exponential Integrals. 


C. Bon 
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Review 44[X, Z], MTAC, v. 11, 1957, 1. 10, p. 48. 
for A. H. Tabu read A. H. Taub. 


‘ 








